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SYSTEMS FOR PREDICTION, RAPID DETECTION, WARNING, PREVENTION OR 
CONTROL OF CHANGES IN ACTIVITY STATES IN THE BRAIN 

BACKGROUND OF THE INVENTION 

1 FIELD OF THE INVENTION 

The present invention relates to the field of neuroscience for analyzing signals representative of a 
subject's brain activity including signals indicative or predictive of epileptic seizures. More particularly, the invention 
concerns the automated analysis of brain activity signals to detect an activity state and transitions between slates, and 
to delect precursors predictive of a change in the subject's activity state to a different stale 

The invention is based on ideas and research in the fields of mathematics, neurology, statistics and 
engineering which enable the real-time analysis of biologic signals such as the eleciro-enccphalogram (EEG) or 
electro-corticogram (ECoG), by the simultaneous execution of multiple methods. In the preferred embodiment, these 
signals are rapidly, accurately, and automatically analyzed in order to 

1 ) Delect and signal the occurrence of an epileptic seizure in real time for contemporaneously with 
the arrival of the signal at the processor/device), 

2) To predict behavioral changes typically associated with seizures, 

3) To predict seizures by detecting precursors to the onset of the electrographic or clinical components 
of a seizure, 

4) To detect and further analyze epileptiform discharges (spikes), and 

5) To download the detection or prediction outputs to devices lor warning, or therapeutic interventions 
or the storage of data. 

2 DESCRIPTION OF THE PRIOR ART 

Humans and animals have several normal stales of behavior such as wakefulness and sleep, as well 
as multiple sub- states such as attentive wakefulness and REM sleep. Abnormal slates include reversible stales such 
as seizures and progressive states such as dementia 

Epilepsy, a disabling disease, affects 1-2% of the American and industrialized world's population, 
and up 10 10% of people in under-developed countries. Electroencephalography is the single most important ancillary 
test in the investigation of ihis disease. EEG's are recorded continuously for hours lo days in an increasing number 
of cases with unclear diagnosis or poor response lo adequate medical treatment. The amount of EEG data for analysis 
is extremely large (e.g.. 64 channels of data at 240 Hz gives 1 .3 billion data points/24 hr or 2.6 Gigabytes/day) and 
consists of complex waveforms with infinite variations 

Visual analysis of these signals remains (exclusive of this invention) the "gold standard" but it is 
impracticable for continuous EEG interpretation as this is the most time -consuming part of any elcctrodia gnostic test 
and requires special training and skills which make this procedure expensive and thus of limited access and use 
Valuable EEG data is often discarded unexamined. The length of recording is unnecessarily prolonged in a specially 
equipped hospital suite until patients have several seizures. If the patient is unaware of the seizures, a common 
occurrence, then a nurse or relative must observe and document the occurrence of these changes. As seizures are brief 
and previously considered unpredictable, the need for continuous observation becomes imperative, adding to the cost 
in a non-effective manner. 
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Present methods of seizure detection are not only expensive, but rely on poorly discriminating 
methods, increasing the review time and nursing assistance because of the large number of false positives, and increas- 
ing the length of hospitalization through the false negatives. Funhermore, these methods often "detect" the seizure 
well after its onset or even its end, when prevention or abatement of the seizure is not possible or irrelevant. 

The inability to process data in real time has thwarted the scientific and clinical development of the 
fields of epilepsy and electroencephalography. Cardiology has developed into a clinical science largelv based on the 
power of electrocardiography to analyze the heart's electrical activity in a rapid and accurate manner This has resulted 
in pace-makers, implanted defibrilators. and other devices which have saved thousands of individuals from premature 
death. The comparison between cardiology/EKG and epilepsy/EEG must take into account the fact that the electrical 
brain signals are far more complex than those originating from the heart. This explains in large pan the developmental 
lag between these two disciplines. 

Electrical brain signals, because of their spatial and temporal characteristics such as non-station- 
arity, have resisted accurate real-time automatic manipulation. The prior an methods presently used lo characterize 
these states are severely limited. For example, the prior an consists of a long history of failed attempts 10 identify 
1 5 changes in EEG during certain behavioral states or tasks and lo discern epi -phenomenology from phenomenology, 
a distinction that would help answer quesuons of fundamental importance. Other limiiaiions include the inability to 
determine whether spikes arc a static marker of epilepsy, or whether they are dynamically related to seizure generation. 

Present methods of automatic EEG analysis have many major limitations which render them 
virtually useless for widespread, safe and effective clinical applications. These limitations include: 
20 1 ) Lack of speed. The time it takes most methods to analyze the input signals and produce an output 

which detects or predicts a state change is too great for use in warning, intervention, or prevention of epileptic seizures 
and other abnormal brain states. 

2) Lack of accuracy. Prior art methods have a large number of false positives (incorrectly iaentifying 
non-seizure activity as a seizure) and false negatives (failure to identify a true seizure), increasing the technical and 

25 financial burden. 

3 ) Lack of adaptability to subject or seizure type, no compromise between speed vs. accuracy. 

4 ) Lack of portability and implant ability 

5) High cost. 

Accurate and reproducible prediction of behavioral or biologic signal changes associated with 
30 seizures has not been possible as these events occur unpredictably. Our methods and devices enable seizure prediction 
by providing a worthwhile prediction time that makes warning, abortion/abatement, and prevention of seizures 
possible. The new treatment modalities that can be based on this method will lead to a significant reduction in seizure 
frequency and consequently, to a reduction in the occurrence of mjuncs and fatalities, allowing persons with epilepsy 
to become productive and lead normal lives. 
35 The prior art in automated seizure and spike detection consists of variations of two primary 

methods: "rule-based" analysis and, more recently, analysis by artificial neural networks. The most popular is a "rule- 
baaed" method which has been under development since the laic 1 970's, primarily by Dr Jean Gotman In the Goiman 
method, the signal is initially replaced by a piecewise linear approximation which connects maxima and minima. 

In the Gotman method, there is a list o! rules which are ihcn applied lo throw out some of the 
40 smaller line segments in an attempt to remove fast activity that is superimposed on an underlying wave of interest. 
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The larger line segments which remain are called "half waves." Gotman's algorithm then compares properties of the 
half waves such as averages of amplitude, duration, rhythmicity, and sharpness in moving 1/3 sec. windows to those 
of past and future data segments As currently implemented, the method uses a total of 30 seconds of past data and 
8- 10 seconds of future data in these comparisons. A set of rules and thresholds are given to determine when these 
comparisons of past, present, and future properties yield a detection of a spike or seizure. 

These rule-based methods have a number of limitations, including a large false positive rate, and 
usually a long delay to detect even abrupt changes (often 10 or more seconds). 

Another method for spike and seizure detection involves training an artificial neural network (ANN) 
using past data tracings with annotated spikes and seizures to "learn" to recognize similar changes in unseen data The 
large number of "neurons" required for accurate analysis of a multichannel EEG/ECoG input signal precludes real -lime 
analysis. Consequently, the current state of the art implementations rely on a smaller number of "neurons" and a 
parametrized input signal in place of the raw signal. The Gotman half- wave decomposition mentioned above is 
commonly used in this signal parametrization step - causing the inclusion of many of the limitations inherent in this 
method to adversely affect the ANN methods In addition, the adaptation of an ANN to improve its performance for 
a particular individual or group is performed off-line and requires time consuming retraining by experienced 
epileptologists This important limitation is overcome by the present invention 
3. GLOSSARY OF TERMS AND USEFUL DEFINITIONS 

The onset of the clinical component of a seizure is the earlier of either ( I ) the time at which the 
subject is aware that a seizure is beginning (the "aura"), or (2) the time at which an observer recognizes a significant 
physical or behavioral change typical of a seizure. 

The onset of the electrographic component of a seizure is defined by the appearance of a class of 
signal changes recognized by electroenccphalographers as characterisiic oi a seizure. This analysts requires visual 
review of signal tracings of varying duration, both before and after the perceived signal changes, using multiple 
channels of information and clinical correlates. The precise determination of the onset is subject to personal 
interpretauon. and may vary based on the skill and attention level of the reviewer, the quality of data and its display. 

The electroencephalogram, or EEG, refers to voltage potentials recorded from the scalp. EEG will 
encompass any recordings outside the dura mater. The elecirocomcogram, or ECoG, refers to voltage potentials 
recorded intracranially, e.g., directiy from the cortex EKG is the abbreviation for electrocardiogram, EMG for 
clcctromyogram (electrical muscle activity), and EOG for elecirooculogram (eye movements). 

The period of time during which a seizure is occurring is called the iclal period. (Those skilled in 
the an will appreciate that the term ictal can be applied to phenomena other than seizures. ) The penod of time when 
the patient is not in (he state of seizure, or in transition into or out of the seizure state, is known as the mien eta I period. 
The preictai period corresponds to the lime of transition between the intcncial and the beginning of the ictal period, 
and the postictal period corresponds to the time period between the end of the seizure and the beginning of the 
interictal period. 

Herein the term real-time describes a system with negligible latency between input and output. 

The term false positive refers to the case of a system mistakenly delecting a non-seizure signal and 
classifying it as a seizure. The term false negative describes the case in which a true seizure goes undetected by a 
system. Systems that have a low rate of false positive detections are called specific, while those with a low rate of false 
negative detections are called sensitive. 
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The terms epileptiform discharge and spike are used interchangeably herein to refer to a class of 
sharply contoured waveforms, usually of relatively large power, and with duration rarely exceeding 200 msec. These 
spikes can form complexes with slow waves, and can occur in singlets or m multiples 

The terms epileptologist and elect roencep ha lographer are used interchangeably. 
SUMMARY OF THE INVENTION 

The preseni invention solves the problems and overcomes the limitations of prior an. while provid- 
ing pioneering advances in the stale of the an lis preferred embodiment enables ( 1 ) the accurate, automated, real-lime 
detection of seizures, as well as the determination of ihcir site of origin, propagation path and speed through regions 
of the brain, and their duration and intensity; (2) the prediction of the onset of the clinical component of seizures; (3) 
the prediction of the onset of the elecirographic component of seizures; (4) the online self-adapiauon. or offline 
adaptation of (1-3) to each individual. (5) the automated use of (1 -3) for diagnosis, quantitative analysis, imaging, 
warning, treatment, and storing of data; and (6) the miniaturization of the system lo a portable or implantable device. 

The adaptation of the system to each individual lakes into account, seizure type and location, and 
changes m the stgnal(s) over time, making use of any existing preictal, ictal. or postictal "fingerprints' for the subject. 
The speed of analysis and levels of sensitivity and specificity can also be adjusted to desired levels, and the method 
can be implemented in either digital or analog form (or a combination). 

The preferred embodiment of the invention uses intracranial or scalp electrodes to obtain signals 
representative of current brain activity and a signal processor such as a personal computer or micro -processor, for 
continuous monitoring and analysis of these signals, detecting imponant changes or the appearance ol precursors 
predictive of an impending change, as soon as they occur. The output of this analysis is then fed to a device which 
produces an immediate response (e.g.. warning, treatment or storage) to the change or predicted change in state The 
signal processing includes an adaptive analysis of frequency, energy, wave shape and dynamics, phase relationships, 
measures of rhythmicity, "scquency.'' and lemporo-spaual siereorypia, variability, dimension, complexity of the signal, 
and noise reduction. 

I . METHODS FOR REAL-TIME SEIZURE DETECTION: 

The following is an overview of the steps which comprise the preferred embodiment of the 
invention for real-time seizure detection. 

(1 ) Extract from the entire signal the pan with ictal characteristics. This step is accomplished with 
adaptive filtering, allowing the selection of patient- and/or sensor-specific initial parameters, and an adaptation process 
in which these filters automatically improve as imponant signal characteristics are learned online. 

(2) The output of this filter is used to compuie an index of ictal activity in each current signal epoch 
(the foreground"), which is then divided by a corresponding measure associated with the background signal, forming 
a rauo. Novel application of median filtering and time and state weighted averaging arc used in this step. 

(3) When the value of this ratio reaches a particular threshold level, a seizure detection is immediately 

signaled. 

(4) Grading and verification of seizures is then accomplished using an analysis of duration, intensity, 
pattern recognition of spatio-temporal propagation, and postictal seizure signal changes. 

In addition, a new seizure imaging method has been developed based on the detection methodology 

presented here. 
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2. METHODS FOR DETECTING PRECURSORS TO SEIZURE: 

These embodiments delect the occurrence of signal characiensiics or patterns which may be 
precursors to the clinical and/or elecirographic components of seizure, resulting in their prediction Determination 
of the onset of a seizure by visual analysis (which is considered "the gold standard") is a subjective and empiric 
process, at the present level of scientific development. Determination of time of seizure onset depends in pan upon 
the specifications and parameters associated with the recording devices and of the location and tvpc of sensors in 
reference to the ussue from where the seizure originates The intensity and degree of spread of the seizure also affects 
detection. 

From a practical standpoint, prediction based on the detection of seizure precursors or the 
elecirographic component itself yield a worthwhile time during which warning and intervention can be instituted to 
abon or prevent the onset of either of the components of the seizure. By virtue of their self- tuning ability (detailed in 
a later section), the continued application of these prediction methods to a given individual or seizure type may 
improve the reliability of subsequent predictions, and may lengthen the prcdictively worthwhile time 

Seizure precursors include, but are not limited to: 

( 1 ) certain significant patterns of epileptiform discharges (or spikes), 

(2) significant abrupt attenuation of signal energy on some or all sensors, and 

(3) significant changes in various characteristics of the power spectral density associated with each of 
the signals that arc being monitored, e.g., a sudden significant drop in the median frequency of a given signal. 

Prediction of seizures may occur during different stages of their temporal evolution: 

a) Prediction of the vibratory (or first) state, i.e., the state before the seizure spreads beyond the 
anatomical or functional boundaries of the "critical epileptogenic mass" (defined as the smallest mass that, once fully 
synchronized, consistently generates the next states). 

b) Prediction of the elecirographic component of seizure. This component is mainly defined by 
temporal continuity of the ictal signal with or without evolution across a frequency spectrum and with some degree 
of propagauon outside the critical mass. Prediction of this state can be made by identifying precursors (see examples). 
Precursors have temporal, spectral, and other characteristics which distinguish them from the elecirographic 
component. 

c) Prediction of the clinical component of seizure. The real-time detection of the electrographic 
seizure component is akin, for partial or secondarily generalized seizures to the prediction of the clinical onset as there 
is a latency between the two components. Precursor detections further lengthen the predictive time of the clinical 
component. 

3. A METHOD FOR SPIKE DETECTION, CLASSIFICATION, AND COUNTING: 

The invention also includes a new method for measuring signal "sharpness," which we refer to as 
least-squares acceleration filtering (LSA- filtering), thai is used as the basis for a new spike detection method h can 
also be used to improve existing spike detection methods. In this method, the signal is continuously monitored for the 
occurrence of such things as spikes, and their amplitude, frequency, waveform, "sequency" (degree or pattern of 
clustering), rale of occurrence, location, and polarity, are computed and stored. This information is then analyzed for 
conformance to any seizure precursor pattern. 
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4. DEVICES FOR THE DETECTION OF SEIZURES. PRECURSORS, AND SPIKES AND FOR THE 

PREDICTION OF SEIZURES: 

The algorithms listed above and defined in detail herein can be realized in digital or analog form 
in a signal processor. The preferred embodiment has been implemented in an Intel 486 based PC for real-time 
5 monitoring and storage of data for patients undergoing clinical evaluation. 
The real-time detection of: 

(a) seizure precursors and the resulting prediction of the electrographic and cluneal seizure 
components, 

(b) the electrographic component and the resulting prediction of the clinical component, or 

10 < c > spikes, enables the institution of safety and therapeutic measures, and initiates or continues the 

adaptation and self-learning of the methods. For example, a seizure prediction can be used to trigger a device for 
systemic, intraventricular, or intracerebral administration of a medicament or substance, for electrical, magnetic, or 
thermal activation or deactivation of a nerve or a region of the subject's brain, for activation or deactivauon of 
physiologic receptors, for ablation of a region of the subject's brain, for activation of a warning or biofeedback device, 

15 or for selection of segments of signals for transmission or storage (or for annotation of continuously recorded signals) 
and further off-line analysis. 

BRIEF DESCRIPTION OF THE DRAWING FIGURES 

Fig. 1 is a schematic illustration of the preferred apparatus of the preseni invention showing inputs 
of brain (or other biologic system) signals of a subject from surface and/or implanted (e.g., intracranial) sensors to a 
20 signal processor and various types of outputs. 

Fig. 2 is a graphical illustration of the segments of data which may be used in one preferred seizure 
detection method for operating the apparatus of Fig 1 to represent current ('foreground') signal aciiviry (e.g., the mosi 
recent 2 seconds), and signal background activity (e.g., a segment of 20 or more seconds in length) delayed 1 second 
from the end of the foreground window; 
25 F, 8 3 A is a graphical illustration showing a finite impulse response (FIR) filter which may be used 

in the first step of the preferred embodiment for detecung seizures in the input signals to the apparatus of Fig. 1 . 

Fig. 3B is a graphical illustration showing the power spectral density (PSD) associated with the FIR 

filler ofFig.3A; 

Fig 4 is a graphical illustration of an ECoG signal used as an input to the apparatus of Fig. I , 
30 Fi 8- 5A is a graphical illustration of the result of applying the generic FrR filter of Fig. 3 A to the 

signal of Fig. 4, and squaring the output signal values m an embodiment of the invention; 

Fig. 5B is a graphical illustration of the dimensionless ratio of a I second foreground median filter 
and a 20 second background delayed median filter applied to the squared, FIR-filiered signal shown in Fig. 5 A in an 
embodiment of an invention; 

35 Fi g- 6A is a graphical illustration of the pan of the ECoG signal from Fig. 4 during which the 

clinical and electrographic seizure onset and subject activation of the event button occurred, and during which the 
apparatus of Fig. 1 detected the seizure; 

Fig. 6B is a graphical illustration of the output of the seizure detection embodiment as presented 
in Fig. 5B but restricting the time window to that time period corresponding to the signal of Fig. 6A. 
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Fig. 7 is a graphical illustration of an embodiment for a seizure imaging method and apparatus 
applied to 16 simultaneous ECoG recordings containing a seizure event; 

Tig 8A is a graphical illustration of an ECoG signal which contains a number of epileptiform 

discharges, 

Fig. 8B is a graphical illustration of the output of first step of the spike detection embodiment; 
potential spikes with sharpness as computed by LSA filtering) exceeding a given threshold and their respective 
polarity (up or down) are identified. 

Fig 9 is a graphical illustration showing the power spectral density (PSD) or the impulse response 
of the presently preferred IIR filler designed to be used m detection of the seizure precursor spikes as described in 
Example 1 

Fig. 10 is a graphical illustration of an ECoG signal which contains a seizure precursor signal (the 
quasi-periodic epileptiform discharge precursor detailed in Example 1 ). In this Figure, the times of the electro graphic 
onset of the seizure, the clinical onset of the seizure, and the lime of detection of the seizure precursor using an 
embodiment of the invention are annotated; 

Fig. 1 1 A is a graphical illustration of an ECoG signal of a subject which contains a seizure 
precursor signal (the signal attenuation precursor detailed in Example 2), 

Fig. 1 1 B is a graphical illustration of the output of the precursor detection embodiment lor this type 
of seizure precursor . 

Fig. 12A shows the illustration of Fig. 1 1 A with the time axis restricted to a smaller range of times 
which still include the detection time of the seizure precursor, and the limes of the electrographic and clinical onsets 
for the seizure; 

Fig. » 23 shows the illustration of Fig. 1 1 B with the time axis restricted to 2 smaller range of times 
which still include the detection time of the seizure precursor, and the times of the electrographic and clinical onsets 
for the seizure; 

Fig. 1 3 A is a graphical illustration of an ECoG signal of a subject which contains another seizure 
precursor signal (a rapid drop in median frequency, as detailed in Example 3). The times of the electrographic onset, 
the clinical onset, and the event button press by the subject are annotated. 

Fig. 1 3B is a graphical illustration of the result of applying the embodiment of the invention for 
detection of this particular seizure precursor to the signal presented in Fig. 1 3 A. The electrographic and clinical 
seizure onsets, and the time of event button press are also annotated, as well as the time of precursor detection; 

Fig. I4A is a graphical illustration of 4.27 sec. of ECoG signal data. 

Fig. 14B is a graphical illustration of the power spectral density (PSD) of the signal in Fig. 14 A, 
showing the modal frequency, median frequency, and mean frequency of the signal; 

Fig. 15A is a graphical illustration of 2 seconds of ECoG signal recorded from a subject during an 
1 nten eta! (non-seizure) period; and 

Fig. 1 SB is a graphical illustration of the absolute value of wavelet coefficients from levels 1-4 
obtained by applying the fast wavelet transform to the signal of Fig. 1 4 A. 
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

Fig. 1 illustrates the preferred apparatus 10 for receiving and analyzing signals representative of 
a subject's brain activity and for producing different types of outputs. Apparatus 10 includes signal processor 12, 
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inpuis 14 and outputs 16. Signal processor 12 is preferably a compuier such as one with capabilities thai meet or 
exceed those of an Intel 4 86 -based computer having 33 MHz clockspeed and 8 MB of RAM. Those skilled m the an 
will appreciate that an appropriate digital signal processor can be used in place of the preferred computer, as could 
a custom-designed semi-conductor chip having the requisite capability, preferably configured for implantation or as ' 
a portable device. Signal processor 12 could also be an analog processor or an analog/digital combination. Appendix 
1 is incorporated as part of the disclosure hereof and contains the presently preferred computer program for use by 
apparatus 10 and, in particular, signal processor 12. and for implementing the preferred methods of the present 
invention as further described herein. 

Inputs 1 4 include EEG (or other type of scalp) signals obtained from a plurality of scalp sensors 
1 8 transmitted through associated lines 22, or ECoG signals obtained from implanted sensors 23 and transmitted 
through associated lines 24. The input signals used in the development of this invention were amplified and convened 
from analog to digital form at a rate of 240 Hz with a dynamic range of [-300,300] ^ V and digital resolution of .59 
(10 bits of precision per datum). This provides 144 Kb of data per minute, per channel Those skilled in the an 
will appreciate that sampling may be performed at fixed or varying rates (higher or lower than 240 Hz) and precision 
(with more or less precision than 10 bits), using linear or nonlinear analog to digital conversion, and with constant 
or varying dynamic range (i.e., adjustable gain) Data acquisition may also be performed using adaptive sampling 
techniques, in which these sampling parameters vary over time and are determined by characteristics of the signal 
being sampled. Adaptive sampling techniques can be used to selectively enhance relevant signal characteristics and 
increase signal quality and resolution in certain frequency bands. 

Outputs 16 can trigger portable or implanted devices, electrodes 26 which may be intracranial or 
extracranial, or placed over or around a nerve 28, a medicament injector or pump 32. an audio or LED output or any 
other form of warning 34, and auxiliary memory 36 for storing input signals and event data Implanted electrodes 26 
can be used for any form of activation or deactivation (e.g.. electrical, thermal, etc.) of local or remoie brain cells or 
for ablation of the epileptogenic tissue Nerve stimulator 28 is preferably associated with the vagus nerve as such 
stimulation has been found to abate or prevent a seizure. Physiologic (or natural) stimulation to receptors (e.g., light 
to retinal receptors) can prevent or abate seizures and is ihe function of stimulator 30. Injector 32 is preferably 
implanted for automated instantaneous release of the appropriate medicament (inclusive of any efficacious substance) 
for treating, preventing or abating a seizure. Memory 36 is provided to store signal and event data for archival and 
analysis purposes. 

As discussed further herein, the analysis performed in signal processor 1 2 can be customized for 
a particular patient to improve the detection of brain states and state transitions, and the prediction of changes in brain 
states. The customization of the signal processing can be based on the information stored in memory 36 via feedback 
of this informauon to signal processor 1 2. For example, this information may be used to monitor efficacy of treatment 
and lo optimize seizure/spike detection and prediction, and therapeutic or safety interventions. Those skilled in the 
art will also appreciate that memory 36 can be included as an iniegral pan of signal processor 12. 

Those skilled in the an will recognize that changes in cerebral state are highly correlated with 
changes m level and type of activity of other organ systems (e.g., heart, etc.) and as such these signals, may be useful 
for detection and prediction or validation of seizures or of other changes in brain slate. The following signals (not 
annotated in Fig. I ) may be used in conjunction with EEG and ECoG signals to further improve performance of the 
system: 
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1 ) Non-electrical cerebral (global or regional) signals, such as concentrations of glucose, free radicals, 
metabolic by-products, neuro-iransrnitiers, or other substances, or measurements of intracranial pressure, temperature, 
blood How or indices of metabolic activity, etc , 

2) Cardiovascular signals such as heart rate, R-R interval and variability, eic, 

3) Respiratory signals such as tidal volume, peak-to-peak interval, etc.. 

4) Elecirodermal and other DC potentials. 

5 ) Signals representative of concentrations in the blood or other peripheral tissues of gases, substances, 
or chemicals such as lactic acid, etc., 

6) Signals representative of the level or type of activity of cranial or peripheral nerves (e.g. frequency 
and partem of action potentials, etc.), 

7) Signals related to EMG activity, force, direction, and patterns of limb or body movements. 
REAL TIME SEIZURE DETECTION 

Successful real-time detection of seizures depends on the ability of any method to rapidly and 
accurately distinguish the ictal from the non-ictal pan of the signal. We begin with a detailed description of the generic 
method, then discuss additional features and modifications which are used to enhance its sensitivity and specificity by 
making it adaptive, i.e., allowing it to learn online. The preferred embodiment as detailed here is based on a sampling 
rate of 240 Hz with 10 bits of precision However, there is a wide range of digitization techniques which may be used, 
together with the appropriate modifications to the algorithm s parameters. Fig. 4 shows a 4 minute segment of ECoG 
signal which is used to illustrate a preferred embodiment for detecting the electrographic component of a seizure. 

The first step of the method consists of applying a filler to the signal to decompose the entire signal 
into its ictal and non- ictal components. As a result of research performed for this invention, identification of key 
differences between ictal and non-ieia! signal characteristics was successfully accomplished These results enabled 
the construction of filters to accomplish this first step of the method. These include "generic" digital filters of both 
finite impulse response (FIR) (also known as a convolution filter and moving average (MA) filler) or infinite impulse 
response (OR) One such FIR filter is shown in Fig. 3, together with an estimate of its power spectral density (PSD). 
These fillers could include analog filters as well and were constructed using a daia base of over 100 seizures in the 
following manner: 

1 ) Each seizure was divided into segments according to its temporal evolution, and the PSD of each 
segment was computed; 

2) The resulting PSD's were then compared with PSD's obtained from lnterictal segments. Power- 
frequenci' envelopes were then computed, more heavily weighing the spectra at frequencies which yielded the greatest 
scparauon between the ictal and interictal PSD values 

3) The orders and rype (e.g. FIR or HR) of the "generic" fillers were then chosen, taking into account 
the trade off between computational burden/speed and goodness -of- fit of their impulse response to the desired shape. 
The filter coefficients were then computed from the desired impulse response using standard filler design methodol- 
ogy 

Those skilled in the an are aware of the many degrees of freedom or options available in designing 
filters. For instance, the magnitude and phase specifications of the filler's impulse response may be adjusted to match 
the PSD characteristics of a sei of seizures, or infinite impulse response (instead of FIR) litters may be used when 
speed of processing is the primary objecuve, especially if the filter must precisely discriminate between extremely low 
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frequencies (e.g.. less than 2 Hz.)- Genetic algorithms provide a useful approach to solving the mulii -dimensional 
constrained optimization problem of designing a best filter for a given set of data 

For a given subject, a filter is initially selected from the filter bank. This selection can be done off- 
line by a clinician upon a study of archival wave forms. In the alternative, the selection can be completed on-imc by 
applying all the filters in the filter bank individually to the input signals and selecting that filter that produces the 
greatest differentiation 

If the input signal is given by fx,;;.,, and the FIR filter has m coefficients { D(h b, b^} t then the 

output (filtered) signal { Y,};., is obtained from the formula 

Yko-twt-b.x^ + b,*^ + -+b 1 ..,x fc ^ M . 
where h is assumed that x, = 0 for all j < I . 

The output of this filter is then squared, and the evolution of the resulting values is monitored, 
comparing the most recent values (the -foreground") to less recent values (the -background"). 10 detect relevant 
changes Fig. SA shows the graph of the Y? values which result from the application of the FIR filter described above 
to the signal of Fig. 4, and then squaring the output. 

In the next step of the method, we shall refer to the present information as "foreground,* and 
compare it to the signal history or -background" Once the raw signal has been filtered to extract and enhance its ictal 
part, the next step in the method is lo use the { YJ} sequence to create measures of the level of ictal acu viry in the most 
recent signal epoch, and compare it with that contained in the interictal signal. To compute a measure of ictal activiry 
in the foreground, we apply a (nonlinear) median filter of order 480 (2 seconds) to the { Y k '} sequence. 10 obtain the 
sequence {F k } (240 values per second per channel). 

F t -median{Y; V(1 Y-^ 3 Y^.Yi*. 

where p is the order of the median filler (e.g., p = 480 @ 240 Hz.). This step is used to monitor a change in central 
tendency of the distribution of the recent (foreground) Y\ values. The median (or other order-statistic is preferred for 
this step because of its ability to measure the central tendency of a distribution without being overly sensitive to 
outliers such as those which result from single or multiple spikes. 

To compute a measure of the ictal activity in the background, (B k which is- used as a reference 
against which foreground changes are measured, we apply another median filter to the Yf values, with the following 
modifications: 

1 The order is increased (e.g.. to 20 seconds, p=4800) to obtain a more stable background, 

2. The filter is delayed by a certain time (e.g., I second) from the current time to remove any possible 
effect of the most recent signal on the background value. B fc , thus making a foreground change easier and faster to 
distinguish/detect, 

3. Updates of background, are disabled during seizures or other anomalies (e.g., transient artifacts). 

and 

4 An exponenualry forgetting time average (or a more general time- and slate-weighted average) of 

this delayed median filter output is used to increase the period of signal history represented by the background 
sequence. (BJ, without increasing memory/storage requirements or computations. Use of this technique decreases 
the size and number of fluctuations in this background sequence, improving detection by allowing more sensitive 
thresholds to be set without significantly increasing false detections. Details regarding more general time- and slate- 
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weighted averaging techniques which may be empioyed in this step are presented in Appendix 2 which is incorporated 
as part of the disclosure hereof 

To now be more precise, the background sequence. (B k } , is computed as follows We begin by 
applying a delayed median filler of order h and delay d 10 the squared output of the FIR-fikcr, { YJ) , to obtain the 
output sequence, { w k } , 

w k = median* YJ^.Y}^.* YJ^.Y^}. 

with, e.g., h = 4800 (to use 20 seconds of data), and d - 240+480 - 720 (i.e., this median filler is delayed 1 second 
from the end of the 2 sec. foreground window, i.e., 3 seconds from the current time), so that 

w k = mcdian{ YJ.jj, Yj.„,.Yj. Tao }. 

Then define 

**' ^ \ B k ifr k *C 2 

where A. is a "forgetting factor" which controls the rale of exponential forgetting (where this factor is between 0 and 
1 ), r k is the current ratio of the foreground "ictal index" to the background "ictal index." given by 




and C, is a threshold level used to disable updates of B k when C, is exceeded by r k . For accurate, real-time detection 
of seizures, the ictal signal should not be allowed to affect the background signal from which the seizures are 
distinguished Accordingly, to further improve the detection method, the constant, C 2 , above was introduced to disable 
updates of the background ictal index, B^, whenever the ratio, r k , is large (i.e., exceeding C : ). Failure to disable these 
updates would allow unwanted rapid growth in the background icial index during a seizure, resulting in a loss of ability 
to measure seizure duration, and a degradation in the ability to detect latter seizures (or the regeneration of a single 
evolving seizure) when these events are closely spaced in time This enhancement of the method is especially 
important when one wishes to detect the end (and hence the duration) of a seizure in addition to the onset time. 

A seizure detection is immediately signaled when the ratio, r k , exceeds C„ which is the detection 
threshold. When the foreground and background median filters are 1-2 sec, and 20 sec, respectively, nominal 
preferred values for the above defined constants are C, = 20, C, - 5, and X ~ .9997. Also note that because the ratio, 
r^ is dirnensionless, the threshold may be set without regard for the units used in the particular signal being monitored. 

Figs. 5B and 6B show the graph of this ratio, when the method is applied to the signal of Figs. 4 and 6A. 
In this example, the method detects the seizure 5 seconds prior to the electrographic onset (as determined 
independently by a trained epileptologist), 8 seconds prior to the clinical onset of the seizure (as determined by review 
of the videotape record by the same epileptologist), and 1 1 seconds prior to the patient's activation of an event button 
(signaling the beginning of the seizure according to the patient), as Fig. 6B illustrates. 

By varying the lengths of background and foreground, the accuracy and speed of detections may 
be adjusted to improve performance as desired for a particular application of the method. The preferred embodiment 
consuiuies a substantial improvement in the state of the art in both speed and accuracy of detection. In addition, the 
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adaptability of the system s parameters provide flexibility and improved performance for a variety of applications. For 
example, if the speed, but not the accuracy of detection is the overriding factor for success in seizure abortion, then 
a decrease in the foreground window length is one way to achieve this result. If, on the other hand, this method is to 
be used in a device to destroy the part of the brain from where the seizures originate, then maximal accuracy (and not 
detection speed) is needed to avoid damaging healthy tissue. This can be accomplished by increasing the length of 
the foreground window. 

The novel use of order-statistic filtering in this invention and, specifically, a median filter, to 
characterize the ictal activity in the foreground significantly improved the accuracy of seizure detection by enabling 
disenrninarion between ictally organized and non-organized epileptiform activity. The further additions of time- and 
state -weigh ted averaging, and the formation of a ratio comparing the level of icta] activity in the recent signal 
(foreground) to that which is normally present during interictal (non -seizure /background) periods, among other ideas 
_ (described below) have enabled the present invention to provide improved results not only in speed, but also in 
accuracy and adaptability. Additionally, the method allows for online learning and adaptation as discussed in later 
sections. 

* 5 As -the signal from each individual sensor is monitored in the aforementioned manner for the 

detection of seizures, spatio-temporal correlations between the ratios simultaneously computed for multiple sensors 
which may be obtaining signals of the same class (e.g., all ECoG), or different classes (e.g.. ECoG, EKG. and blood 
cheniistry) can then be used, if necessary, to increase sensitivity and specificity of the algorithm, e.g.. by allowing for 
sensor and signal dependent threshold settings and the use of mu Hi -channel information in both seizure detection and 

20 artifact rejection. Those skilled in the an recognize that external sensors produce a lower signal to notse ratio than 
implanted sensors. This drawback is particularly prominent for sensors placed on the scalp for EEG monitoring due 
to volume conduction and low-pass filter effects on the conical signal caused by the structures surrounding the brain. 
Furthermore, the signal recorded from the scalp may be delayed from that directly recorded from the brain. In order 
to compensate for these drawbacks, the following strategies have been adopted (individually, or in combination): 

25 1 ) Prefiliering of the scalp signal and other inputs to extract the useful signals (e.g., separate conical 

yollage potentials from artifacts). 

2) Artifact pattern recognition. Artifacts, defined as signals that are not the direct product of cortical 
neuronal activity occur during normal activity and also during seizures. During their clinical component, seizures 
generate a host of artifact signals that, for a given individual, may be highly stereotypical in their spectral and spatio- 

30 temporal characteristics and as such are distinguishable from cortical seizure activity. These artifacts, which 
correspond to body, mouth, or eye movements, etc., are first recognized by comparison with artifacts that have been 
catalogued into a general library according to their characteristics. As these artifacts are identified, their temporal 
evolution is tested for conformance to a pattern of artifacts stereotypic of that individual s seizures. 

3) Use of other signals such as EKG. respiratory rate, magnetic fields, thermic potentials. 
35 concentrations of chemicals, recordings from a dynamometer or acceleromeier attached to the patient or bed, etc. for 

use in seizure detection or prediction. 

Once a seizure has ended, there is usually a decrease in signal power in certain frequency bands. 
This decrease is dependent upon the intensity and duration of the seizure, among other factors. The method/device 
tests for and measures the duration of any loss of power following a large increase (such as that which occurs with 
40 seizures), and the results of these tests are used to retrospectively assess the accuracy of seizure detection. For 
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instance, a large sustained increase in power, without a subsequent power decrease in certain bands, would be 
routinely subjected to off-line review. 

THE ADAPTIVE AND EVOLUTIONARY NATURE OF THE METHOD 

Large scale validation studies have proven that the above-detailed embodiment tor seizure detection 
is both highly sensitive and specific. However, to account for possible intra- or inter-individual signal variability, and 
to further improve performance as needed, a number of additional steps have been implemented which aiiou the 
system to learn and adapt on and offline. These steps can be grouped into the following categories; 

( 1 ) adaptive signal acquisition methods. 

(2) intelligent parameter adjustment, and 

(3) adaptive use of detection and prediction information. 

The first step in adapting the method to a particular subject or application consists of selecting the 
appropriate set of signals to be monitored, and controlling the manner in which these signals are acquired Although 
the detailed example above makes use of signals consisting of electrical potentials (EEG/ECoG) sampled at 240Hz 
and digitized with 10 bits of precision, in some cases one may, as discussed earlier, vary the analog to digital 
conversion parameters, or use the analog signal itself, in order to ascertain various characteristics of the signal which 
may be important in the subsequent analysis. For example, sampling rate may vary continuously with the frequency 
content of the signal, increasing with the slope of the wave (differential sampling); the steeper the slope, the higher 
the sampling rate. In addition, online signal quality control methods can be used to detect various forms of signal 
degradation and warn the user or others. For example, the system can produce an output indicating the presence of 
a large quantity of 60 Hz activity, or abnormally high "clipping" in the analog to digital conversions, etc Moreover, 
in many cases it is advantageous to record signals from other sites (or organs) in addition to the brain or its encasing 
structures, and to analyze biological signals other than EEG/ECoG, including, but not limited to, ( 1 ) other electrical 
potentials (AC or DC), (2) EKG, (3) EMG, (4) respiratory rate, (4) concentrations of glucose, free radicals, or other 
substances (e.g.. neurotransmitiers) in the brain, blood or other tissues, (5) magnetic fields, and (6) brain temperature 
It is also important to note that while, for the sake of clarity, attention is primarily focused on data from a single sensor 
in the detailed example above, in practice the method can be applied to signals from a large number of individual 
sensors and combinations of sensors (e.g., a weighted sum of the electrical potentials between a given sensor and all 
other sensors) simultaneously (i.e., in parallel), monitoring spatio-temporal correlations of the resulting outputs. 

The parameters of the method may be adjusted as needed Parameters such as foreground and back- 
ground window length, filter shape and type, time- weight, and threshold settings may be adapted on or offline to the 
particular subject and application. For example, if data exists containing a prior seizure for a given subject, then one 
can process this data to realize a filter adapted to any known signal characteristics or seizure "fingerprint" of that 
subject. A filter with the spectral characteristics of its impulse response matching the typical PSD ai seizure onset 
for this subject, can be used to initialize the adaptive filtering procedure, thereby improving the sensitivity and 
specificity of the algorithm for that subject The FIR filter may also be adapted to the particular location of the sensor, 
or the type of signal being monitored. For example, when monitoring posterior electrodes (e.g occipital), the 
preferred filter is designed to recognize and de-emphasize irrelevant signals in the alpha range (8- 1 3 Hz with normal 
reactivity), when its power is below a given percentile While spectral characteristics of seizures for a given subject 
may differ from those analyzed for the design of the generic filler, this method, by virtue of its adaptability, is well 
suited to prediction and detection of seizure patterns with a wide range of spectral and other relevant characteristics 
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To increase compuLauoaal efficiency-, one may also use a stable infiruic impulse response (IIR) filter 
( also known as an auio- regressive moving-average or *ARMA" filler) in place of the FrR filter. Such fillers use a 
linear combination of past filtered signal values in addition to present and past input signal values to compute each 
new filtered signal value. That is, given IIR filter coefficients 

5 A=[a, a,aja, ... aj, and B=[b, b, bj... bj, and input signal {x,, x 2 x K } . we compute the sequence { Y N ) using the 

recursive formula: 

This feedback of the filtered signal allows the IIR -filtered output to be produced with fewer computations and shorter 
delay time. 

_ The FIR filter step utilized in the preferred embodiment example is a special case of more genera) 

1 0 adaptive filtering methods which may be employed to decompose the signal, extracting and enhancing signal activity 
mdicauve of an ongoing or impending change of brain state. For example, one may allow the FIR filter coefficients 
(and order) mentioned above to vary with time as a function of the current brain state, to increase the accuracy of the 
method under certain conditions. The method can use a different filter when the subject is in slow- wave sleep, from 
that which is used during vigorous exercise, since the characteristics of the background signal and noise level/type, 
1 5 from which a change is detected, arc markedly different in these two stales. Those skilled in the art are aware of many 
known methods of digital signal processing which can be used to achieve such online adaptation (e.g.. "active noise 
cancellation") In the presently preferred implementation, one begins with a generic filter (as above), that allows for 
the detection of seizures which may begin in a wide range of frequency bands. As the input signal evolves over time 
in the non-seizure state, windowed power spectra] density (PSD) estimates of the signal can be successively computed. 
20 and a PSD representative of the recent (or entire) signal history may be obtained, together with confidence intervals, 
a median, min, max, and other percentiles, etc. This representative PSD is then used to modify the current filter's 
impulse response, in accordance with the newly acquired subject -specific and state-representative "background" 
. information to improve detection of subsequent state changes, lniericial PSD's which do noufit any predetermined 
ictal panem, but which are sufficiently different from the background PSD, may be archived and reviewed; these and 
25 parameters which may be computed from the PSD (see Appendix 2), may be used as templates for delecting other 
state changes, precursors to stale transitions, and for signal quality control. Archived seizure segments can also be 
used in the online adaptation of this filter, focusing attention on the frequency bands of past seizures which are 
maximally differentiated from their respective lmerictal segments. Those skilled in the an will further recognize that 
several methods exist for online filter design from a given PSD. e.g., the method of Parks- McClelland. In addition 
30 to filters whose design is mainly based on the PDS, phase, shape and other char act ensues (eg, neural lv -based filters 
or other signal -shape fillers) can be used to realize new filters as necessary. 

Parameters can be adjusted in this method, to detect ictal activity on a wide range of time scales and 
sensitivities, in particular, by reducing the length of the foreground window from the 2 second length used above to 
.2 seconds, individual spikes can be detected with this method The detection of these spikes as ihey occur enables 
3 5 the creation of a biofeedback system in which the subjects are made immediately aware of epileptiform activity in their 
brain (eg., by an auditory transducuon of the abnormal signal activity for auditory pattern recognition). Such a system 
can be used to teach the subjects to control or eliminate their abnormal own discharges. 
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Thc thresholds. C, and C,. can also be adjusted online, based or. accumulated siaiisiics and 
functional of ihe compuied rauos, {r k } , over long time periods, including ihe rouo values attained aunnp seizures, 
together with ihe maximum and minimum values reached intcrictally 

Finally, the ways in which ihe system uses detection and prediction information can be made 
adaptive as well. For example, a warning alarm could have a feedback mechanism to detect background noise, and 
adjust its output signal accordingly - from a vibration when in the library to a loud lone when on an airport runway 
As another example, changes in stimulation parameters for intervention or seizure abortion can be made based on a 
history of therapeutic effectiveness. 

OTHER EMBODIMENTS OF THE SEIZURE DETECTION METHOD 

The task of seizure detection requires near- simultaneous temporal -frequency localization of the 
signal. Real-time seizure detection imposes the additional demand of high processing speed. 

The use of adaptive filtering in the preferred embodiment is ideally suited 10 meet these demands. 
Those skilled in the art will appreciate that a variety of other methods can be used for nearly- simultaneous temporal - 
frequency localization, including the windowed Founer(or Gabor) transform (WFFT). the wavelei transform, the 
Hartley transform, and the Wigner-Ville transform. Any of these methods can be used to extract or enhance the ictal 
part of the signal from the entire signal as in Ihe first step of the adaptive filtering method 

In addition, a number of other transforms can be used to accurately and rapidly extract and enhance 
the ictal portion of the signal. We give two examples: 

1 Windowed correlation integrals: Embed the original (windowed) signal in a higher- 

dimensional space using ihe method of time-delays, a standard technique in nonlinear 
dvnamics and statistics. Then count the number of pairs of points whose separation is less 
than some critical distance (this is called the sample correlation integral). Monitor this 
statistic as a function of time, ictal activity is indicated when the number of such pair falls 
by an order of magnitude or more. Higher order correlations can also be used. 
2. The windowed "kinetic energy," defined as follows: take the first difference of the time 

series represented by the signal, re-embed this derived time series with a suitable lag 
time, then choose a window size (time interval), in each window, and compute for each 
window the sum of the squared lengths of all ihe vectors, and monitor this statistic as a 
funcuonofiime. 

Both methods are robust against changes in the control parameters such as embedding dimension, 
delay time, and window length. In particular, both have been effectively employed wilh short window lengths and 
embedding dimensions as low as 3, enabling real-time monitoring, detection and prediction. Both of these methods 
measure whal might be characterized as a relative dispersion, at a certain scale, of the points on the re-embedded 
trajectory As such, they arc remarkably insensitive to small fluctuations in ihe position of the points, which in turn 
means that these methods arc extremely robust against the contamination of the signal by noise. For example, using 
method 1, only a slight decrease in sensitivity occurs when the signal is contaminated by sufficient Gaussian noise to 
reduce the signal -to- noise ratio to 5dB. 

It should be noted that ihe trarisforms listed in this section are also useful for preprocessing signals 
when hnle is known a priori about the frequency bands of interest or the time scale of changes. In particular, these 
methods can be used as initial screening tools in determining frequency bands in which changes are correlated wilh 
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pamcular changes in brain state Additional background details regarding the Founer and wavelet mcihods may be 
found in Appendix 2. 

Coherence analysis, and higher order spectra and cumulanis. etc., also provide additional important 
information regarding signal frequency changes over time. 

The present ly preferred computer program for seizure detection shown in Appendix I performs 
on-line median filtering, updating the moving foreground and background median filters at 240 Hz. For this task, the 
program makes use of circular doubly-linked lists. For certain applications (e.g., to conserve processing resources 
when a large number of signals are being monitored simultaneously), these computations may be performed less often 
using well known batch sorting algorithms (e.g., Quicksort), to compute the median in moving windows. One may 
also replace the median filter by a similar order-statistic filter, or with some other measure of central tendency (e.g.. 
a -trimmed means) For compulation of the background index described above in the case when a very large 
background window is desirable, one may instead compute the sequence, {B^} using the exponentially forgetting 
average of the median of periodically sampled foreground values For example, the background value can be updated 
once per second by first computing the median of the foreground values that occurred ai each of the last 300 seconds, 
then adding this result (properly weighted), to the previous background value. 

Another preferred embodiment of the above methodology consists of the analog implementation 
of the digital methods described herein. Those skilled in the art will appreciate that even' slep of the method presented 
herein can be implemented in analog This fact is important in the context of miniaturization and the development of 
implantable devices, since analog devices generally require significantly less battery power to operate and, thus, last 
much longer. 

AN APPLICATION OF THE SEIZURE DETECTION METHOD TO IMAGING 

The seizure detection method applied to a signal recorded from a particular sensor produces a 
dimensionJess ratio representing the current level of ictal activity present at that corresponding location. As pan of 
the present invention, a new "seizure imaging" method was developed based on this seizure detection methodology. 
The seizure image is obtained using a contour plot of the brain regions which achieve the same seizure ratio ("equi- 
ictaT), and tracking the temporal and spatial evolution of this ratio; the contours are directly obtained by interpolating 
these ratios over both space and time. 

Fig 7 illustrates the application of this seizure imaging method to a seizure recorded simultaneously 
from two implanted needles each containing 8 contacts Contacts 1 -8 recorded activity from the left, and contacts 9-16 
recorded activity from the right amygdalo-hippocampal regions. The times of seizure detection (solid line, 1 07 sec. 
after an arbitrary initial zero time), elecirographic onset (dashed line, 1 12 sec.), clinical onset (dot-dashed line, 1 15 
sec), and event button press (dotted line, 1 1 8 sec ), are annotated This graph illustrates that the seizure originates 
from contact 1 2, then weakens somewhat, followed by a more widespread resurgence at neighboring contacts, first 
in those more anterior (9-1 1 ), and then later involving the more posterior right temporal contacts (1 3, 1 4). The onset 
of activiry in the left hemisphere begins at 1 56 sec. on contact 4, but then spreads to contacts 1-3, and 5 within 3 sec, 
and to the posterior contacts 6-8 four seconds later. Ictal activity on the left is evident for 1 5 seconds after the 
cessation of right temporal involvement. This particular imaging technique can be very helpful to the clinician in 
spat*o- temporal localization of seizure activity and. in particular, in localization of the site of origin of the seizure and 
characterization of the pattern of spread. 
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It has been found as pan of the present invention, through the application of this seizure imaging 
method, that seizures from the same subject usually have a high degree of spa no- temporal congruency of ratio 
intensity Accordingly, such images and the pathways of ictal propagation which they depict, can be used in the 
analysis of spatial correlations in the outputs of the simultaneously evolving, single -channel detection ratios computed 
as described in the preferred embodiment of the method Specifically, any detected seizure may be compared in terms 
of its spatial evoluuon, propagation speed, and intensity evolution (including, e.g., the trajectory of maxima) absolute 
and relative intensities), with similar information from prior ictal events, and the degree of conformance can be 
quantified and used in assessing whether (a) the detection was a seizure, or (b) seizures originate from one or multiple 
sites These methods may also allow therapeutic intervention at the sites to where the seizure spreads, when treatment 
at the site of on gin is of limited efficacy. 

Those skilled in the art will appreciate that other quantities computed from the input signals (not 
just the above-mentioned ratios) can have their spatio-temporal evolution depicted via similar imaging techniques, 
and that a variety of different existing interpolation techniques may be employed. 
A METHOD FOR DETECTION OF EPILEPTIFORM DISCHARGES OR SPIKES 

The method described herein for detection and classification of spikes, is based on a new method 
developed as part of this invention for computing relative sharpness of the signal as it evolves, and comparing it to 
the sharpness of other waveforms. In prior an methods, the sharpness is computed at each successive signal maximum 
or minimum using the numerical second derivative at that point, or by connecting adjacent signal minima and maxima 
with straight lines and measuring the resulting angles. These approaches, however, are inaccurate because they ignore 
signal characteristics near each peak, and rely too heavily on a small number of digital points approximating the 
underlying continuous signal, using too few points near the peak (the region of primary interest) and too many points 
far away from the peak. The present invention includes a new method for computing this sharpness using what we 
shall call a "ieast -squares acceleration filter* (USA filter), which overcomes these limitations of prior an. This method 
can be used independently, or as an improvement to existing methods, many of which rely heavily on relative signal 
sharpness in their detection of spikes. 

This method consists of first fitting a function to the data surrounding the point at which the 
sharpness is to be computed, and then determining the sharpness of the fitted function. The preferred function for use 
in this fit is a parabola, and the preferred criterion for the fit is the minimization of mean-squared error of fit (least 
squares). The parabola is used because its second derivative can be obtained from the computation of a single value 
(its leading coefficient), and no function evaluations (which can be computationally expensive) are necessary. Other 
performance criteria for an optimal fit could also be used. The choice of least squares is relatively standard, highly 
effective, and has been made computationally feasible here. 

The computation of the sharpness (i.e., the acceleration) of the signal(s) at each data point has been 
reduced to applying a small order FIR filter to the data from that sensor. The output of this filler is equivalent to first 
computing the parabola which attains the optimal fit (in the least squares sense) to the data near the point in question, 
then computing the second derivative of this parabola, which is used as a measure of sharpness. 

One may then obtain an index ratio of relative sharpness at a particular point in the signal by 
dividing this computed peak sharpness by, e.g., the median sharpness of all waves occurring in a moving background 
window of the signal. If the absolute value of the resulting dimensionlcss relative sharpness exceeds a given threshold. 
C„ then the wave is classified as a potential spike. The preferred nominal value for C s is 50. but may be adjusted as 
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needed When a pdeniial spike is detected, the polarity of the discharge may then be obtained from the sign of the 
computed acceleration. 

Fig. 8 shows a 10 second segment of the ECoG signal containing a number of spikes, along with 
the resulting detections (and their respective polarities) made using this method. Now spatio- temporal correlations ' 
5 may be used to (a) reject artifacts, (b) group detected "spikes" into complexes, and (c) detect the occurrence of 
discharge patterns which may serve as precursors to seizure. Source code useful in accordance with the present 
invention is included in the microfiche appendix. Those skilled in the an will appreciate that genetic algorithms alone 
or in conjunction with neural network or fuzzy logic may be effectively used to analyze the spatio-temporal correla- 
tions of detected potential "spikes" for artifact rejection. 
10 PREDICTION OF SEIZURES THROUGH THE DETECTION OF PRECURSORS 

It has been found as part of the present invention that many patients have precursor signals that 
_ consistently precede the onset of both the electrographic and clinical components of their seizures. Three examples 
described below illustrate such precursors, present detailed methods for detecting them in real-lime, and demonstrate 
their utility in seizure prediction. 
15 EXAMPLE 1 

Often subjects have an onset of large energy, primarily low frequency, quasi -periodic epileptiform 
discharges which occur seconds to minuies prior to the electrographic component of a seizure Such patterns arc rare 
to non-existent in the interictal state The following is a description of how the present invention may detect the 
occurrence of such a precursor: 

20 First, a linear combination of the signals from relevant sensors is formed in a manner that enables 

the preservation of the polarity of each discharge so that there is no cancellation through superposition of waves. Then 
a filter designed to extract precursor spikes with patient -individualized shape is applied to this composite signal. A 
generic choice of niter for use in this step (when the signal is sampled at 240 Hz), is an I1R filter with PSD as shown 
in Fig. 9. This filler was designed using a data base of these precursor signals from various patients in a manner 

25 similar to that described earlier for the generic seizure detection filter. An IIR filter is preferred here instead of an FIR 
filler because of its advantages in filtering into such low frequency bands as that required here. The output of this filter 
is then squared, and, as in the seizure detection method, a background signal is computed by applying a 20 second 
median filler to these squared values, and then exponentially forgetting (slowly) the result. The relevant spikes are 
then identified by detecting the instant that the ratio of the current filtered signal (squared) divided by the background 

30 signal value exceeds a threshold value. C«. Since these spikes are now individually detected, the method then tests 
for their occurrence in a particular rrrythmic manner which is highly correlated with a later onsei of the electrographic 
and then dmical components of seizure Specifically, we impose "periodicity constraints" which tesi if there is at least 
t t seconds and at most U seconds between spikes, and lhat at least N spikes occur in succession according to these 
spacings before a detection is signaled. The preferred nonadaptive settings of the above constants are C, = 100, t, = 

35 l,u= 10. and N = 2. Each of these parameters (and the filter employed) can be adapted to patterns which may known 
a priori, or learned on-line for a particular subject via retrospective analysis of previously detected seizures. 

Fig 10 presents a graph of an ECoG segment recorded from one of the subjects thai exhibits this 
particular precursor, and the time of precursor detection using ihis method is annotated, along with the onset times 
of the electrographic and clinical components of the seizure. This method has produced a prediction of clinical onset 

40 an average of 54 seconds prior to its occurrence, and a prediction of clcctrogrnphic onset an average of 42 seconds 
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prior to its occurrence. The detection of this precursor has been followed by a seizure within two minutes in 1 00% 
of cases analyzed to date. 

example: 

As pan of the present invention, it has been found thai signal attenuation or "quieting" precedes, 
5 by up to tens of seconds, the onset of the clinical and electrographic components of seizure in many subjects In those 
that exhibit this preictal quieting, the application of the following preferred 

method to detect the onset of this attenuation results in a prediction of the electrographic and clinical seizure 
components. 

For any signal, {X, , i > 0}, the average signal energy over the interval of time, t, < t < u, is given 

10 by 

A{t v t 7 y— x — p x; ds. 

'i J/ . 

The average signal energy, E„ in a moving lime window ot length T seconds is given by 

E, - A(,-T.,) - i f ' T X* ds. 

which can be computed recursively and efficiently online using the formula 

£,*a, = £, * A(U+&t) - A{t-T.t + bt-T) 

or 

The long-time average energy, H,, used as an adaptive background value against which energy changes are measured, 
15 is given (recursively) by 

with the preferred value of lambda being slightly less than I The above recursions can be initialized using the 
formulae: 



Now a ratio, R,, is computed as 
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and a precursor detection is made the instant that this ratio exceeds a threshold value C.. Note that an increase in this 
ratio corresponds to a decrease in average signal energy. Preferred nominal values for the parameters in this method 
are T = 5, and C 5 = 5. but these again may be adapted to particular subjects. This method has produced a prediction 
of clinical onset an average of 23 seconds prior to its occurrence. • 

Figs. 1 1 and 12 present a graph of an ECoG segment recorded from a subject that exhibits this par- 
ticular precursor, and indicates when the detection is made using this method relative to the clinical and electrographic 
onset limes of the seizure. The detection is made on the graph 278 seconds after an arbitrary zero, which is 15 
seconds prior to the independently determined electrographic onset, and 1 9 seconds prior to the time of clinical seizure 
onset 

EXAMPLE 3 

It has been found as part of the present invention that, for some subjects, certain sudden changes 
in the power spectral density (PSD) of the signal may be used to predict an impending seizure onset. For example, 
a sudden significant drop in the median frequency of the signal (defined in Appendix 2) ts a consistent precursor to 
seizure for some subjects. The following is a description of the preferred method for detecting the occurrence of such 
a precursor: 

Begin by computing the median frequency of the particular signal of interest in moving windows 
of length T, (as described in Appendix 2). Compute a background median frequency using an average of the median 
frequency values in a moving window of length T ? . Then compute the ratio of the background median frequency 
divided by the median frequency in the current T, second window. If this ratio exceeds a threshold value C , a 
detection of this precursor is immediately signaled. The preferred nominal/non- adaptive choices ol* parameters are 
T, = 2567240 sec. (approx. 1 .067 sec ), T, = 60 sec, and C 6 = 5. 

Fig. 1 3A shows a 5 minute segment of ECoG data which contains a seizure in the last minuie. The 
times of electrographic and clinical seizure onset, and the lime at which the event button was pressed are annotated 
(dashed, dash-dot, and dolled lines, respectively). Fig. 1 3B shows the graph of ihe ratio described above over this 5 
minuie segment, with a detection occurring 1 2 seconds prior to the electrographic onset, 1 5 seconds pnor to the 
clinical onset, and 1 8 seconds prior to the patient's pressing of the evem button. In the example used in Fig. 1 3, it 
should be noted that there is some signal energy attenuation as well (Ihe precursor described in the previous example), 
but that this is not used in the detection of the precursor described now 

Other precursors which may be relevant for predicting an impending seizure for a particular subject 
or group can be isolated using modem pattern recognition techniques. As mentioned previously, preictal or lnterictal 
patterns which are present in the output of our seizure imaging method (which uses spatio-temporal interpolation of 
the ratios or other relevant quantities computed at each sensor) arc strong candidates for precursors 

Another useful method for determining precursor signals for a given subject or group makes use 
of waveform classification and cluster analysis, followed by patiem recognition techniques applied to the resulting 
sequence of waveform "Labels." The following illustrates how this can be done. One begins by segmenting an arriving 
input signal into individual "waves.* This can be done in a variety of ways, for example, one could define a wave as 
5 1 2 cansecuuvc data points, or, more adoptively, as the number of data points between 2 successive baseline crossings 
of ihe signal. As each new waveform is segmented, it can be extensively analyzed and subsequently classified into 
a tpossiblv evolving) library of distinct waveforms or "cyma.** 
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c 

For example, the duration, amplitude, number of base line crossings, arclenglh, maximal sharpness, 
number of local maxima and minima, area bounded between the waveform and a horizontal axis, total power, and the 
fraction of power in a particular frequency band (and any number of other quantities) could each be computed (and. 
if desired, statistically normalized), with the resulting n measures used to form a vector in n -dimensional space 
5 quantifying that particular waveform (in this example, n=9). The distances (either Euclidean or non- Euclidean) 
between this new "point" in n -dimensional space and a library list of other "points" can then be computed and analyzed 
to see which waveform or "cymcme" in the library most closely resembles the new waveform. In this way, the input 
signal can be labeled as a sequence of cymemcs (or indices into a list of waveforms), and the resulting "cyma" can be 
analyzed as described above for the occurrence of patterns which may serve as precursors to a given change of state. 

1 0 Such waveform lists can be constructed from available data using currently known methods of cluster analysis, and 
neural networks can be employed in making online or offline decisions classifying a given waveform Specifically , 
we have successfully used a competitive learning network trained with the Kohonen learning rule and an adaptive 

- learning rate which decreases over time for this task of classifying the sequence of n -dimensional points into groups. 
These analyses can be performed consecutively where, in each step, a different analysis method can be used 

i 5 CORRELATION DETECTION AND AUTOMATIC PRECURSOR IDENTIFICATION AND ISOLATION 

This section describes the preferred methods for automatic identification of signal characteristics 
which arc significantly correlated with later changes of state in this subject's brain. This is based on correlational 
analysis and pattern recognition techniques applied to the signal recorded prior to each change of state. The 
transformations which apply to a given input signal or set of signals (e.g., FIR filtering, fast wavelet and Fourier 

20 transforms, etc.) decompose the input and present its information content in an organized and useable form. The 
original signal, and the product of these transformations are then automatically analyzed for the occurrence of patterns 
significantly correlated with detected changes of state 

Signal analysis for correlation may occur on-line, or off-line from signals previously recorded and 
stored in memory Further, the analysis can be done on the basis of predetermined patterns or patterns developed 

25 through correlation analysis, or both. One may select the segments of signal for analysis which precede these changes 
of state with "markers" in some externally controlled way, or one may simply let the software accumulate correlations 
between "major changes" and precursors in an inclusive fashion for eventual use. It is also important to note that the 
process can either be done off-line using data obtained from the subject or putative patterns, or on-line by automated 
systems installed as part of the device. 

30 The following is an example describing the use of correlations for determining average signal 

power. Let b ft be the I th time coefficient in the k* wavelet basis expansion (see next section for more details), using 
a fundamental wavelet time scale dt- Let the average value of b tt over an interval, chosen lor illustration to be one 
minute, be denoted < b t (l) >. This is an average over roughly 1^= 2 U * values of 1. The interval begins one minute 
prior tt> time t. Let < p(t) > be the average power in the signal over the same time interval, and < p > be the average 

35 power over many previous intervals, such as 1 00. Form the running deviation of coefficients from the average, b tt - 
<b,(t)> 

Form the running deviation of the power from its long term average. < pf t) > - < p >. A correlation 
matrix C k (k) is formed by multiplying 

C h (k) = (<p(t)> - < p >Xb a -<b k (t)>). 
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This matrix depends on the level, k. the umc step within the interval. 1. and the interval start time t If there is no 
correlation between the signal and the power in a particular level k, then the matrix for thai k value will be a random 
variable in the index I. If a persistent correlation begins at umc i*+ l*dt, at some level k*. then the product will rise 
at point 1* and stay large, up to statistical fluctuations- 
Applying this procedure over many intervals of "normal" or background state, typical values of 
C„(k) will be fairly independent of 1 and t. and only depend on the level k. Distributions of these values are formed 
and evaluated When C>(k) rises wcU above a slausucally likely value it represents a correlation between b u and the 
power, < p(0 >. in that interval. In one application, restricting the analysis first to intervals marked t* where < p(t*) 
> rises weU above the background, one identifies the region of the seizure. Only the C b (k) values for the immediately 
preceding interval need to be kept in this case. Then by examining the correlations C„.(k) as a function of level k and 
time 1. the level or levels k* where large correlations occur for I > 1* can be determined automatically. The lime for 
the precursor to predict m advance the seizure is found by computing (1^ - |*)dt. These values of k» and I* can be 
output and stored, leading to automatic recognition of precursors 

Another similar method is to iram any assembly consisting of a number of logic elements with 
connected threshold rules (such as a neural network) to recognize correlations between the precursor signals, and the 
macroscopic change of state, in this case the correlations between precursors and state change may be automatically 
registered and analyzed without the intermediate step of forming correlation matrices 

The same technique can be applied to linear combinations of the signal after application of various 
signal filters, e.g.. considering simultaneously derived wavelet coefficients from different levels of resolution after 
application of the discrete wavelet transform. One could also project the input signal(s) into various other bases (i.e., 
apply other signal decompositions) and monitor the resulting signal components as they evolve over time ("pattern 
recognition"). 

THE USE OF GENETIC ALGORITHMS AND GENETIC PROGRAMMING IN THE ADAPTATION AND 
EVOLUTION OF THE METHODS FOR SEIZURE. SPIKE, AND PRECURSOR DETECTION 

Algorithms (GA) and genetic programming (GP) may be used as part of the overall strategy of 
algorithm adaptation and evolution in the methods for seizure and spike detection, and seizure prediction. The 
problem of detection and prediction of seizures is a highly complex one, due to the large intra- and inter- individual 
variability of brain signals imenctally and ictally. The problem of selecting parameters such as filter coefficients, 
thresholds, window lengths, type of signals, etc. to optimize a particular performance criteria is ideally suited for the 
use of GA In this method, parameter sets are assembled into chromosome- 1 ike strings which are interchangeable and 
form the "first generation." This first generation is evaluated by a fitness function which will assign a score to each 
dement of the generation. After the scoring is completed, the "best" elements are saved and reproduced for the next 
generation, possibly combined together ("cross-over"), and/or new elements ("mutations") are introduced. This 
process of generation, reproduction and testing is repealed, until elements of the resulting generation achieve the 
desired level of performance for the particular application. 

The following example illustrates how GA may be utilized to choose a time-in variant FIR that 
maximizes performance. The first generation consists of a group of bandpass filters, each with 2 Hz bandwidth, of 
varying orders (but less than 240 Hzl), centered at odd integral frequencies between 1 Hz. and 1 1 9 Hz., and designed 
according to the Parks- McClelland method. The resulting sensitivity and specificity arc measured when each filler 
is used in the seizure detection method presented herein (e.g.. by computing the mean squared time between detection 
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and ihc electrographic onset as indcpcndcnlly scored by an electroenccphalographer), and the best fillers are saved 
and "reproduced" for the next generation. These fillers are also combined together via cascading ("cross-over"), and 
new "miuation" filters of different band widths are added at random frequencies. Mutations may also mvoive randomly 
decreasing or increasing the order of certain filters from one generation to the next, but staying below the upper limit 
of 240. and/or using another filter design methodology. This evolution is continued until some stabilization is reached, 
or until the resulting performance achieves e high enough degree of sensitivity and specificity for the given data base. 

Gencuc programming, a form of program induction, enables the method/system to evolve (without 
external inputs or re-programming) based on internal fitness functions GP may enable the method to self- learn and 
extract the most relevant and useful characteristics of the signal. 

Having thus described the preferred embodiments of the present invention, the following is claimed 
as new and desired to be secured by Letters Patent: 
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APPENDIX 2 

SOME BACKGROUND INFORMATION ON MATHEMATICAL METHODS 

This Appendix presents some background detail helpful in understanding the preferred methods 
of the present invention. 

THE DISCRETE FOURIER TRANSFORM 

It is well known thai every discrete signal defined on evenly spaced discrete time points { t„. 

, J . where l> = 2nk/N. can be interpolated by a trigonometric polynomial which can be written as 



where 



' 6 =0j\/ = (A' - l)/2 ifN is odd 
, 6 = !>/ = (N - 2)/2 // N is even 



and 



. AM 
N k=0 



This discrete Founer transform (DFT) represents a given signal (or "time series") as a superposition 
1 5 of sine and cosine waves with different frequencies and amplitudes. The number of Founer coefficients is equal lo 
the number of data points in the signal and reproduces (in some sense) the information contained in the original signal. 
In practice, the DFT is computed using the fast Fourier transform (FFT). Once the signal is decomposed into these 
fundamental frequencies, one may analyze the signal's components in separate frequency ranges (bands), and examine 
the coefficients for dominant frequencies. One generally computes an estimate of the power spectral density (PSD) 
20 from the Fourier coefficients in order to determine the dominate modes of the system The simplest such estimator 
is the periodogram. which is a graph of the squares of magnitude of each of the Fourier coefficients, in order to sec 
dominant modes of the system. As used herein, PSD means any of the estimators commonly employed in signal 
analysis. 

The PSD of the signal {f(to),f(l,) t ...fO^,)} is obtained from the Founer coefficients, {c,J. as 



\c/r if j = o 

2!c| : i/0 < j < M ♦ 6 - I 
'/J = A/ ♦ 9 •- I . 
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Hcre p, mav be interpreted as ihe total power at frequency w ( in ihc segment of signal which was trans: ormcd There 
are M f 0 different frequencies. w r at which signal power is computed, and these frequencies are evenly spaced 
between w 0 = 0 H2. and w MMM = F/2 (the so-called Nyquist frequency ) which is half of the samplmc frequency, F,, 
of the signal (240 Hz above). The PSD contains precise frequency information about the given signal ("frequency 
localization"), but contains no inform at ton about the particular times al which a given frequency occurred However, 
in the preferred applications, the time of a particular frequency change is important 

The windowed fast Founer transform (WFFT) was created to address the lack of temporal 
resolution of the discrete Fourier transform, the signal is partitioned into successive segments or windows and the PSD 
in each window is computed, allowing the user to track changes in the frequency content of the signal window-by- 
window, i.e., in time. It must be understood, however, that the temporal resolution of the lime of frequency changes 
in the signal is only on the order of the length of the window. Also, due to the discrete nature of the data, there is a 
degradation in the ability to accurately compute the full frequency content of a signal as the size of window is 
decreased since the number of frequency coefficients returned by the DFT is equal to the number of data points 
transformed (the window length). 

An important application of Fourier methods in the present invention involves the analysis of the 
temporal evolution of various parameters associated with the frequency domain representation of the signal (i.e., the 
Fourier transformed signal). To accomplish this, one first computes the windowed fast Fourier transform (WFFT) of 
the signal obtained from a single sensor in a moving window of current data values and then the corresponding PSD 
estimate for each window Then each PSD curve is converted to a probability density function (pdf) via normal izauon, 
and the resulting densities are treated as if they were characterizations of a non-stationary random frequency. To be 
more precise, suppose p(x) is a non-negative function defined on the interval, [0,L],(0 < x < L) Without loss of 
generality , one may assume that the area under the graph of p(\* \? eq::s! ;c one. : .e . one rr.cy trs.r r :r - >± r !f 
fjpp(x)d\ * 1 , then one may simply redefine p via normalization, as 

p(x) 
Jo p{x)dx 

Considering p(x) as a pdf, one can compute statistics of the signal frequency distribution using 
probabilistic methods. In doing so, we obtain a number of parameters from the pdf in each window of time and 
monitor their absolute and relative temporal changes along with the time evolution of various relations between them. 
Some of the parameters computed are measures of central tendency in the PSD, including order statistics such as the 
median frequency (w.r.l. PSD), various moments of the pdf (e g , the mean frequency), modal frequencies (i.e., the 
frequency with maximum power); and other parameters measure fluctuations in frequency in each window (e.g., ihe 
frequency variance, skewncss, and kurtosis) 

Precise formulae are now presented for some of these measures which arc most commonly used. 
If p(x) is a continuous p.d.f. (defined discretely in this case) which is zero for x < 0 and x > 1.. then the mean of the 
distribution (also known as the "first moment." or "expected value") is given by 
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L 

u = f xp(x)dx. 
o 



the XI th moment is given by 

L 

0 

The variance of the distribution is given by V- The median, nv of the distribution is defined by the equation 

"V 

| p(x)</r - 1 - 50%. 
o * 

Note thai one may also define other percentiles (order statistics) in this manner. The median frequency of a signal in 
a given window of data is the frequency which divides the area under the PSD in half The mode, M f , of the p d.f. is 
5 given by 

M r = argmax{p(x)}. 

The modal frequency of the PSD (i.e., the frequency at which the maximum power occurs) is thus defined. Fig 1 4 A 
presents a signal comprised of 1024 points of ECoG data (about 4.267 sec.) recorded from a subject during an 
intenctal (non-seizure) period Fig. 14B illustrates the corresponding power spectral density (PSD) of the signal, 
showing the modal frequency, median frequency, and mean frequency of the signal, computed according to the above 
definitions. 

One can use the modal frequency and variants of this concept to test the signal for rhythnuciry, i.e., 
to detect segments of data which are nearly periodic. When this quasi -periodicity of the signal occurs, the power in 
the signal is concentrated at a few resonant modal frequencies. Many useful measures capable of detecting 
hypers>Tichronous patterns of neuronal firing often found in the recruitment and entrainment associated with seizures 
can be derived using a combination of measures such as those defined above. 

For example, to detect hypersynchronous behavior of neurons focusing attention near 1 5 Hz activity, 
one can compute a frequency-biased functional of the modal frequency such as 

where p is the PSD of the signal in, for example, a moving window of 256 points, and M f is the modal frequency for 
that window. One may then monitor the evolution of this measure for significant increases relative to background to 
produce a rapid detection of this type of signal acti viry 
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This method for deiccung rhythmic discharges in ihc signal can be further enhanced by weighing 
a measure such as thai above by other measures of quasi- periodicity thai can be computed from the signal itself 
(without first applying the FFT). For example, the reciprocal of the standard deviation often successive lmer-zero- 
crossings (or interpeak intervals) increases dramatically in the event of hypersynchrony. The product of this measure 
with the one mentioned in the last paragraph provides an excellent method for the detection of this type of phenomena. 

A number of other nonlinear functions based on these quantities (e.g., the product of the average 
energy and the inverse of the median frequency) can also be utilized as a means to obtain precursor information 
regarding an imminent change of brain state. 
THE DISCRETE WAVELET TRANSFORM 

Whiie the Fourier transform described above gives precise frequency localization of a particular 
signal, the discrete wavelet transform (DWT) has recently gained popularity because of us ability to provide useful 
temporal and frequency information simultaneously The DWT expresses a given signal in terms of a basis consisting 
of translations and dilations of a single "mother wavelet," W(j). To be more precise, a given signal is 
expressed as 

1 IJc 

here 

-[ 

The wavelet coefficients are defined by 

/V 

The first "dilation" index of the coefficient controls the level of resolution of information, while the 
second "translation" index controls the temporal information contained in the coefficient. Representing a given signal 
in terms of wavelet coefficients rather than in the usual lime domain is analogous in many ways to representing a 
musical composition as sheet music rather than looking at the graph of the sound waves through an oscilloscope as 
the piece is played. Each chord on the sheet music contains information about both the frequencies that are played 
and the time at which they occur. 

As with the DFT, fast, efficient algorithms exist to compute the wavelet coefficients { b a J . The fast 
wavelet transform (FWT) which is based on the pyramid algorithm developed by Mallat makes the FWT even easier 
to compute than the FFT. 

The FWT is applied to a given signal in windows of data of length 2" points, for some prescribed 
value of n (e.g., n=6). The window size may vary based on a particular application and implementation Use of the 
FWT may be desirable if one wishes to monitor signal changes in multiple frequency bands simultaneously and these 
desired frequency bands match those obtained by dilations of a particular mother wavelet. One may then group the 
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wavelei coefficients by level," i.e., by like dilation factor, ordered in each level according 10 their progression in lime. 
The lime between level 1 wavelet coefficients is 2 l dt. where di is ihe sample interval of the original signal Hence, the 
level 1 coefficients contain the finest resolution (i.e.. the highest frequency) information contained in the signal, and 
each lower level provides signal information at increasing time scales and decreasing frequency ranges. Level 1 corre- 
sponds to the smallest scale of resolution, and contains twice as many coefficients as level 2, which contains twice as 
many as level three, and so on. Recall thai the EEC signal used in illustrating the preferred seizure detection 
embodiment described above was sampled at 240 Hr, so that a new data point appears roughly even' .004 seconds. 
Computing the FWT using, e.g.. 64 data points (.267 sec.) results in 32 level I coefficients, 1 6 level 2 coefficients, 
8 level 3 coefficients, etc. Thus for each level, there is a corresponding lime series of wavelet coefficients representing 
the temporal evolution of signal power in a number of frequency bands, each higher frequency level "covering" a 
frequency band of width twice that of the next lower frequency level Fig. 1 5 shows (the absolute value of) these series 
_ for levels I -4 for a typical inlenctal segment of 512 data points (about 2 seconds) together with the original signal. 
A LEAST SQUARES ACCELERATION FILTER 

In the method described earlier herein for computing the sharpness of a given waveform at a 
1 5 particular point, one first obtains an optimal parabolic fit to the data near the point in question One then may use the 
second derivative (i.e.. the accderauon) of the resulting parabola as a measure of sharpness at the peak. The criterion 
of minimizing the mean square error of this fit may be applied lo obtain the best filling parabola tin a least squares 
sense). This subsection provides the necessary mathematical details for accomplishing this step. 

^U>v>\). j - 1 nj are data points to he interpolated by the parabola. p(x) = a,\ : + a,x + a^, ther. 

the optimal coefficients a 2 , a„ and ao which minimize the mean square error 



10 



20 



£ (v - /><* )>-. 
7 = 1 



over all choices of such coefficients, are obtained by solving the equation 

K a 0 y = cv.\V(.v*n. 



where X, = xf. for 1 < 1 < n. 1 < j < 3. Y = [y,...yj\ and X 4 denotes the Moorc-Pcnrase pscudoinverse of the matrix 
X. If we let A = X*X and B = X*Y, then we have A^ = x^\ 1 £ i, j * 3.and the column vector B, = xJX 
1=1,2,3. Once one solves this cquauon to obtain the optimal value lor the leading coefficient, a ;i the second derivative 
of the parabola is equal to twice this optimal coefficient The sign of this second derivative indicates the polarity of 
the potemial spike as well, e.g.. a negative second derivative indicates that the peak of the spike is a local maximum 
for the signal. The above formulae for computing a, may be simplified using the symbolic pseudo- inverse of the 
3x3 matnx and subsequent matrix multiplication, resulting in the following formulae for a "least squares acceleration 
filter" which have been implemented to compute this sharpness measure in real time (i.e.. the compulation of the 
second derivative at a given point requires only p(p- 1 ) floating point operations): 
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a(p.di) * di 4 p (p- 1 ) (2p- 1 ) (3p : -3p- 1 )/30, 
b(p.dt)*=di J p=(p-l) J /4, 
ctp.di) = di I p(p-l)(2p-l)/6, 
d(p,dl) = dlp(p-l)/2, 
e(p) = p. 

A * c e - d\ 
B = cd-be, 
C = bd-c\ 

D = aA + 2bcd-eb : -c\ 
E = |ABC]/D, 

f=[Odl : (2dl) : (3dt) 7 ... ((p-))dl) 3 ], 
g = (Odi 2dt ... (p-i)dil, 
h = [ I I \ ... 1], 

F-2E[f g«hT. 

and 

a,(k-U(p- 1 )/2])) = F(p) x(k) + F(p- 1 ) x(k- I ) + ... 

F(2) x(k-p+2) + F(l)x(k-p+l), 
where {x L ,k= 1 ,2,... } is the signal being analyzed, p is ihe number of points used in the parabolic fit (e.g. p=7), and 
dt is the time step of the signal being analyzed. Note thai, for a fixed value of p, the filler coefficients F may be 
computed once and stored for later use in the compulation of a, from the above FIR filter The delay in computing a, 
at a given point, which is instilled by using [[(p-l)/2]] (i.e., the greatest integer less than or equal to (p-1 )/2) future 
data points, is only (p+1 )*dt/2 seconds. For example, with p=7 and dt= 1/240, the delay is 1/60 sec. 
TIME- WEIGHTED AVERAGING 

In the methods described above, there are many cases in which it is desirable to compute some 
background or reference value for a particular signal. By accurately representing the history of the signal, one may 
improve the method's ability to identify relevant changes which standout from this background. 

In this invention, time-weighted averaging is preferred. A subsei of these techniques are able lo 
determine a suitable long time-average of any desired functional of the input signal in a very computationally efficient 
manner, taking into account the entire recorded history of the signal and using only a minimal amount of computer 
rnemory. A particular example of the more general method is that of exponential forgetting, in which the recent signal 
information is weighted more heavity that the distant past. 

The general form referred to as a time-weighted average of a continuous-time signal (X,.t > 0 } with 
time-weight {f 0 ,i > 0, 0 < s < t} is given by {m,,t>0}, where 

,U r . 

' /„* 

0 
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Thc discrete version of this lime-average is obtained by simply replacing the integrals in the above definition by the 
corresponding summations over the index variable s For ccnain limc-weights. the above formula mav be written 
recursively, in particular, this may be achieved for the case when f is independent of t. If the time- weight f u = e", then 
a version of this time-weighted average can be simplified to the exponentially forgetting method that are employed 

5 in some of the embodiments of the invention described herein. 

Variants on this choice may be useful for certain applications. For example, by making X a periodic 
function of j with a period of one day, the time-weighted average may be used to weight signal information at a 
particular time of day more heavily than at other times. This may be particularly important for an individual who may 
commonly have seizure events only during certain times of day. 

0 Also note that the choice of time- weight f tl = X\u+\ 8' vcs to the usual moving average 




Where % denotes the indicator function. 

In the more general case of "time and state weighted averaging," the weight function, f, may also 
depend upon the signal X itself. This additional generalization incorporates the case in which the historical average 
not only depends on the epoch length over which the signal value is averaged, but also upon the signal value over this 
epoch. This technique may be useful when one desires to allow certain signal characteristics to, bv their existence, 
modify the desired signal background average. 
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Appendix i 

The following is a group of KAT1AB scripts and function 
files, and C subroutines which are used in our signal 
analysis. We have omitted some plotting routines to 
display the data. The code is organized by it's particular 
use and a brief description of each piece of code is 
given. 

online. c 

This c program contains subroutines necessary to implement 
the seizure detection algorithm in real-time. There are 
two primary functions, and utility functions that they call 
and a MATLAB Cmex "gateway" program which calls them. 
SetupCompute initializes arrays and pointers at the 
beginning of monitoring, while Compute is called once for 
each acquired raw data point. The program filters the raw 
signal using a FIR filter, computes the median of filter 
coefficients (squared) in a foreground moving window, and 
compares the ratio of this median to one similarly computed 
from a background window delayed from the end of the 
foreground window. If this ratio is above a preset 
threshold, a seizure detection is made. 



**+★★ online. c 

/* Mark G. Frei 10-10-94 */ 



online. c 



The calling syntax is: 

[ratios, fg_vals, bgjvals) - online5_mex (x , _ 
f ilter_coef f s, 



Input arguments: 
x 

f ilter_coef f s 

FG_LENGTH 

DELAY_LENGTH 

BG_LENGTH 

THRESHOLD 



FG_LENGTH , DELAY_LENGTH , 
BG — LENGTH , THRESHOLD) ; 



signal to be filtered; 

filter coef f itients; 

foreground length; 

delay length; 

background length; 
threshold. 



Output arguments: 
ratios 
signal to the 

f g_vals 
bg vals 



— the ratio of filtered foreground 

filtered background. 
filtered foreground signal 

— filtered background signal 
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#inclucfe <stdlib.h> 

#include <math.h> 

/*— ™ #include <stdio.h>*/ 

^include "mex.h" 

/* #include "compute . hpp" */ 
/* #include- n pdmal6 . hpp" */ 

/* extern Status_Word statusWord; */ 
/* extern float ratio; */ 
/* extern int raw; */ 

/* extern Computation_Parameters compParam; */ 

float raw; /★ Should be global */ 

float fg_med, bg_med, ratio ; 

int DATA_LENGTH, FG_LENGTH, DELAY_LENGTH , BG LENGTH , 

NPTS, FILTER_LENGTH , THRESHOLD; 
double *f ilter_coef f s ; 

/* DEFINITIONS: */ 

#define max (A, B) ((A) > (B) ? (A) : (B) ) 
#define min(A, B) ((A) < (B) ? (A) : (B) ) 

#define FILTER_LENGTH_MAX 600 

#define lambda 0.9997 /* Forgetting factor for background 
*/ 



#define bg_threshold 5.0 /* Level of ratio used to control 
when to update background * / 

float new_data ( FILTER_LENGTH_MAX ] ; 

/* Here are the definitions of the data structures */ 

struct data_node{ 

struct data_ node *next; 

struct index_node * ind_ptr ; 

float data; 

} (*filter_data) [ ] ; 
/* Define pointers to this linked list */ 
struct data_node *begin_f g_data , *end_f g_data , 
*begin_ bg_ data , *end_bg_data ; 

struct data_node *begin_delay_data , *end_delay_data ; 

struct index node { 

struct Tndex_node *prev, *next; 

struct data_node *dat ptr; 

} (*ind_fg)[), (*ind_5g)[]; 
/* Define pointers to this linked list */ 
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struct index_node * i__f g_:nin , *i_fg_med, *i_bg_nin, 
* i bg_ined ; 

struct index_node *ind_place, *ind_fg_new f *ind_bg_new; 

Function Prototypes: 
***★■*★**★★**★******★★****★*/ 

void SetupCompute ( void ) ; 
void Compute ( void ) ; 

struct index_node *remove_node ( struct index_node *iptr) ; 
int insert_after (struct index_node *place, struct 
index_node *new_node) ; 

struct index_node *search( struct index_node *iptr, struct 
index_node *min_j?tr, float val) ; 

/*★★****★******★********★★*****★★******★*★+*★★+**★**★★+*★* 

void SetupCompute ( void ) 

/* This function initializes all the data structures and 

arrays used in Compute * / 

{ 

int i; 

/* Allocate big arrays */ 

f ilter_data=calloc(NPTS, sizeof (struct data_node) ) ; 
ind_f g=calloc (FG_LENGTH, sizeof ( struct index_node) ) ; 
ind_bg=calloc (BG_LENGTH, sizeof (struct index_node) ) ; 

/* Initialize new_data array */ 

for (i=0; i < FILTER_LENGTH ; i++) 
new_data [ i ] =0 . 0 ; 

/* Initialize and connect filter_data, ind_bg, ana ind_fg 
structures */ 

bg_nted«*0 . 0001 ; 
ratio=l . 0 /THRESHOLD ; 

for (i=0; i < NPTS-1; i++) 
{ 

(*f ilter_data) (i) .next = k ( ( *f ilter_data) [ i+l] ) ; 
(*f ilter_data) [i] . data = bg_med; /* Some 
arbitrary value */ 
> 

(*f ilter_data) [NPTS-l] .next = & ( ( *f ilter_data) [ 0 ] ) ; 
(*filter_data) (NPTS-l] .data = bg_med; 

for (i»=0; i < FG_LENGTH-1; i++) 
{ 

(*ind_fg) (i+l) .prev = £. ( (*ind_f g) [ i ] ) ; 
(*ind_fg) [i] .next = & ( ( *ind_f g) [ i+i ] ) ; 
(*ind_fg) [i] .dat_ptr « & ( (*f ilter_data) [NPTS - 1 
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ll) ; 

(*f ilter_data) [NPTS - 1 - ij.ind ptr = 
&( (*ind_fg) [i]> ; 

> 

(*ind_fg) [OJ .prev = k ( (*ind_f g) [FG LENGTH - l])- 
(*ind_fg) [FG_LENGTH - lj.next = &(7*ind fg)[01)' 
(*ind_fg) [FG_LENGTH - lj.dat ptr = ~ 

k ( (*f ilter_data) [NPTS - FG_LENGTHj) ; 

(*filter_data) [NPTS-FG LENGTH] . ind ptr = 

i ( (*ind_fg) [FG_LENGTH - 

for (i=0; i < BG LENGTH- 1; i++) 
{ 

(*ind_bg) [i+i] .prev = M (*ind_bg) [ i ) ) ; 
(*ind_bg) [i) .next - & ( ( *ind_bg) [ i+i ] ) ; 
(*ind_bg) [i] .dat_ptr = k ( (*f ilter_data) [ i ] ) ; 
^ (*filter_data) [i] . ind_ptr = k ( (*ind_bg) [ i ] ) ; 

(*ind_bg) [0] .prev k ( (*ind_bg) [BG LENGTH - 11); 

(*md_bg) [BG_LENGTH - lj.next = S(7*ind bgWO))- 

(*indjag) [BG__LENGTH - l } . dat ptr = ~ 
& ( (*f ilter_data) [BG_LENGTH - 

(*f ilter_data) [BG_LENGTH - l].ind ptr = 
k ( (*ind_bg) [BG_LENGTH - 1]); 

/* Initialize other pointers to structures */ 

begin_fg_data = k ( (*f ilter_data) [NPTS - l]); 
end_fg_data = k ( ( *f ilter_data) [NPTS-FG_LENGTH ] ) ; 
begm_delay_data = k( (*f ilter_data ) [ NPTS— FG_LENGTH - 

end_delay_data » 
k ( (*f ilter_data) ( NPTS-FG_LENGTH— DELAY LENGTH) ) ; 

begin_bg_data - ~ 
k( (*filter_data) [NPTS-FG LENGTH-DELAY LENGTH - 1-1); 

end_bg_data = k ( (*f Ilter_data) ( oj) ; 

i_fg_min = k ( (*ind_fg) [0] ) ; 

g.aed «= & ( (*ind_fg) [ ( FG_LENGTH / 2 ) ] ) ; 
x_bg_min = k ( (*ind_bg) [0] ) ; 
i_bg_med « k { (*ind_bg) [ (BG_LENGTH/2 ) ) ) ; 

} 

/****+*+++*******++*+***+**+***+++++*++++++++++++++++++^++ 

void Compute ( void ) 
{ 

int i, iflag_fg, iflag_bg; 

float fc, intjrem, int_rem_bg, int_med, int add, tmpi , 
d tap; ~~ 



bg_med_tmp ; 
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/+ Begin by computing the inner product of recent data with 
filter */ 



fc= raw*( (float) f ilter_coef f s [ FILTER_LENGTH-1 ] ) ; 

for (i=0; i <= FILTER_LENGTH-2 ; i++) 
{ 

fc= fc+ new_data[i+l]*f iltcr_coef fs[i] ; 
new_data [ i ]=new_data [ i+1 ] ; 

} 

newjata [ FILTER_LENGTH-1 J =raw ; 
fc=fc*f c; 

/* Store oldest background filter_data value for a later 
comparison */ 

int_rem_bg= end_bg_data->data ; 

/* Detach "old" index_nodes corresponding to end_fg_jdata 
and end_ bg_data */ 

/* making a temporary copy of each so no pointers are lost 
*/ 

ind_f g_new=remove_node ( end_f g_data->ind_ptr ) ; 
ind_bg_new=remove_node ( end_bg_data->ind_ptr ) ; 

/* If a node corresponding to a minumum value is removed, 
then move min pointer up one */ 

if (ind_fg_new » i_fg_min) 
i_f g_min=i_f g_min->next ; 

if (ind_bg_new == i_bg_min) 
i__bg_min=i_bg_rain->next ; 

/* If a node corresponding to a median value is removed, 
set a flag to recompute it later */ 

if lag_fg«=0; 

if (ind_fg_new — i_fg_med) 
if lag_ f g«=l ; " 

if lag_bg=0 ; 

if (ind_bg_new ==* i_bg med) 
if lag_bg=l ; 

/+ Search ind_fg (beginning at place where value 
corresponding to begin_f g_data */ 
/* was placed) to find where the filter coefficient 
corresponding to end_bg_data */ 

/* (our new filtered data value) belongs in sorted order. 

*/ 
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tmpl= ( i_f g_med->dar_ptr->data ) - (begin_f g_data->data) ; 
tmpl = tmpl* (f c- ( i_f g_med->dat_ptr->data) ) ; 
if (tmpl > 0) 
{ 

ind_ place=search ( i_fg_roed, i_fg_min, f c) ; 

} 

else 
{ 

ind_place=search ( begin_f g_data->ind_ptr , i_fg_min, 

fc) ; 
} 

/* Insert detached node after place determined by the 
search * / 

insert_af ter ( ind_place, ind_fg_new ); 

/* Interconnect the newly placed index node with the 
corresponding data_node */ 

ind_f g_new->dat_ptr - end_bg_data; 

end_bg_data->ind_ptr *» ind_fg_nev; 

/* Update the minimum pointer (if necessary) . Note that if 
the new filter_data value */ 

/* is less than or equal to previous minimum, then the new 
node will be inserted between*/ 

/* i_fg_min and i__f g_min->prev the latter of which has its 

dat_ptr pointing to the */ 

/* maximum foreground filter_data value 

*/ 

if (fc <= i_f g_min->dat_ptr->data) 
i_f g_ min=i_f g_min->prev ; 

/* Update median pointer (if necessary) */ 

int_rero«end_ f g — data->data; /* filter_data value 
removed from foreground */ 

int_med=*i_f g_med->dat_ptr->data ; /* median from 
previous foreground window */ 

int_add=fc; /* newly added filter_data 

value */ 

if ((int_rem < interned) && (int_add > int_med) ) 

i_fg__med=i f g_med->next ; 
if ((int_rem > Tnt_med) (int_add <= int_med)) 

i_f g_med=i_f g_med->prev; 
if ((int_rem = int_med) || (if lag_f g==l) ) /* Update 
"the hard way" */ 
{ 

ind_place«i_f g_min ; 
for (i=l; i<«=(FG_LENGTH / 2) ; i++) 
ind_place-ind_place->next ; 
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i_f g_med=ind_place ; 



/** Now we do the same thing to find the median for the 
background window **/ 

/* Search ind_bg (beginning at place where value 
corresponding to */ 

/* begin_bg_data was placed) to find where the filter 
coefficient */ 

/* corresponding to end_delay_data belongs in sorted order. 
*/ 

tmpl= ( i_bg_med->dat_ptr->data) - (begin_bg_data->data ) ; 
tmpl = 

tmpl*( (end_delay_data->data)-(i_bg_med->dat_ptr->data) ) ; 
if (tmpl > 0) 
{ 

ind_place=search ( i_bg_med, i_bg_min, 
end_delay_data->data) ; 
} 

else 

{ 

ind_place=search ( begin_bg_data->ind_ptr , i_bg_min, 
end_delay_data->data) ; 
} 

/* Insert detached node after place determined by the 
search */ 

insert_af ter ( ind_place, ind_bg_new) ; 

/* Interconnect the newly placed index node with the 
corresponding data_node */ 

ind_bg_new->dat_ptr = end_delay_data ; 
end_delay_data->ind_ptr = ind_bg_new; 

/* Update the minimum pointer (if necessary) . Note that if 
the new */ 

/* filter_data value is less than or equal to previous 
minimum, then */ 

/* the new node will be inserted between i_fg_min and 
i_fg_min->prev */ 

/* the latter of which has its dat_ptr pointing to the 

maximum */ 

/* foreground filter_data value 

*/ 

if (end_delay_data->data <= i_bg_roin->dat_ptr->data) 
i_bg__min= i__bg_min- >prev ; 

/* Update median pointer (if necessary) */ 
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int_rem=int_rein_bg; /* filter_data value removed from 
background */ 

int_med=i_bg_med->dat_ptr->data ; /* median from 
previous background window * / 

if { if lag_bg==i) 

i 

int_med=int_rem_bg ; 

} 

int_add=end_delay_data->data ; /* newly added 
filter_data value */ 

if ((int_rero < int_med) &£. (int_add > int_med) ) 

i_bg_med=i_bg_med->next ; 
if ((int_rem > interned) &6 (int_add <= int_med)) 

i_bg_roed=i_bg_xned->prev ; 
if ((int_rem == int_xned) || (if lag bg==i) ) /* Update 
"the hard way" */ ~ 

{ 

ind_place=i_bg_min ; 

for (i=l; i<=(BG_LENGTH/2) ; i++) 
{ 

ind_place=ind_place->next ; 

} 

i_bg_med=ind_place; 

} 

/* Now place new filter coefficient (squared) into place in 
data_node * / 

/* vacated by the last value in background which we are 
throwing away */ 

end_bg_data->data= fc; 

/★*★★* Move begin/end data^node pointers forward **★**/ 

begin_f g_data=begin_f g_data->next ; 
begin_bg_data=*begin_bg_data->next ; 
begin_delay_data=begin_delay_data->next ; 
end_f g_data«end_f g_data->next ; 
end_bg_data=end_bg_data->next ; 
end_delay_data«end_delay_data->next ; 



/* Compute ratio from medians and return */ 



f g_med= i_f g_med->dat_ptr->data ; 
bg_med_tmp= i_bg_roed->dat_ptr->data ; 
if (fabs(bg_med_tmp) < le-12) 

ratio=0 . 0 ; 
bg_med=bg_med_tjnp ; 
} else { 

if (ratio >= bg_thr eshold / THRESHOLD ) /* don't 
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update bg_med */ 
{ 

ratio - fg_med / bg_roed ; 
} else { /* update bg_med using exp 

forgetting */ 

bg_med = lambda *bg_med + 
( 3.- lambda ) *bg_raed_tmp ; 

ratio = fg_jned / bg_med; 

} 

ratio = ratio /THRESHOLD ; 

) 

} /* end of compute */ 

/+ + + + + + + + + + + + + +*++ir + + + + it + + + + 4+ +* + + * + * + + + + + + + + + + + * + * + + + + + + + * 

struct index_node *remove_node ( struct index_node *iptr) 
/* This function detaches an index node and connects the */ 
/* adjacent nodes, returning a pointer to the detached node 
*/ 
i 

( iptr->prev) ->next = iptr~>next; 
( iptr->next) ->prev = iptr->prev; 
return (iptr) ; 
} 

+ I 

int insert_after (struct index_node *place, struct 
index_node *new_node) 

/* This function inserts the index node pointed to by 
new_node after the */ 

/* node pointed to by place (i.e., the data value" 
corresponding to */ 

/* new_node is slightly larger than that corresponding to 

place) */ 

{ 

new__node->prev=place ; 
new_node->next=place->next ; 
place->next->prev«new_node; 
place->next=new_node ; 
return 0; 
> 

struct index_node *search( struct index_node *iptr, struct 
index_node *min_ptr, float val) 

/* This function performs a linear search through filter 
coefficients */ 

/* corresponding to elements of an index_node, to find the 
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piece */ 

/* where the new filter coefficient, val, belongs in an 
ordered list */ 

/* The input index_node pointer, iptr, points initially to 
where the search */ 

/* begins and is incremented as the search proceeds 

*/ 

/* The function returns a pointer to the index node after 
which a node * / *~ 

/* corresponding to val should be inserted in the list 

*/ 

{ 

float temp; 

temp=iptr->dat_ptr->data ; 

if (val > temp) /* Need to search larger values 

than initial guess */ 
{ 

while ((val > temp)) 
{ 

iptr=iptr->next; 
temp=iptr->dat_ptr->data ; 
if ( iptr==min_ptr ) 
< 

iptr=iptr->prev; 
return (iptr) ; 

} 

} 

iptr=iptr->prev; /* Back pointer up one node since 
we've just passed it */ 

return ( iptr) ; 
) else 
{ 

while (val <= temp) 
{ 

iptr=*iptr->prev ; 
temp«iptr->dat_ptr->data ; 
if (iptr==min ptr->prev) 
return (iptr) ; 

} 

return (iptr); /* No need to back up in this 

direction of search */ 
} 

} 



static 

#ifdef STDC 

void online ( 
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double ratios[] ( /* Declare output */ 

double fg_vals[], 
double bg_vals[]/ 

double x(] /+ Declare inputs */ 

) 

#else 

online (ratios , fg_vals, bg_vals, x) 

double ratios[], /* Declare output */ 

f g_vals [ ] , bg_vals [ ] ; 
double x[]; /* Declare inputs */ 

#endif 
{ 

int i,jj, temp; 

struct index_node *iptr; 

SetupCompute ( ) ; 

for (jj=l; j j<=DATA_LENGTH; 
i 

raw «= (float) x[jj-l); 
Compute ( ) ; 

ratios ] = (double) ratio; 
fg_vals[ j j-l]= (double) fg_med; 
bg_vals( j j-l]= (double) bg_med; 

} 

return ; 

} 

/* GATEWAY FUNCTION */ 

#ifdef STDC 

void roexFunction ( 
int nlhs, 
Matrix *plhs[], 
int nrhs , 

Matrix *prhs[] 
) 

#else 

mexFunction(nlhs, plhs, nrhs, prhs) 

int nlhs, nrhs; 

Matrix *plhs(], *prhs[); 
#endif 
{ 

double *ratios, *fg_vals, *bg_vals, *x; 

/* Check for proper number of arguments */ 

if (nrhs i~ 6) { 

mexErrMsgTxt ("ONLINE requires six input arguments."); } 
else if ((nlhs != l) && (nlhs != 3)) { 
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mexErrMsgTxt ( "ONLINE requires one or three output 
argument . " ) ; 

> 

F I LTER_ LENGTH = max (mxGetM ( prhs [ 1 ] ) , mxGetN ( prhs [ 1 ] ) ) ; 

DATA_LENGTH = max (mxGetM (prhs ( 0 ] ) , mxGetN (prhs ( 0 ] ) ) ; 

/* Create matrix for the return argument */ 

plhs [ 0 ) « mxCreateFull (DATA_LENGTH , 1 , REAL) ; 
ratios = mxGetPr (plhs [ 0] ) ; 

if (nlhs ==3) { 

plhs[l] = mxCreateFull (DATA_LENGTH , 1, REAL); 
fg_vals = mxGetPr (plhs[ 1 ] ) ; 

plhs [2 ] = mxCreateFull (DAT A_LENGTH , 1, REAL) ; 
bg_vals = mxGetPr (plhs ( 2 ] ) ; 

} 

x = mxGetPr (prhs [ 0] ) ; 
f ilter_coef f s = mxGetPr (prhs [ 1 ] ) ; 
FG_LENGTH *= ( in t ) ( mxGet Pr ( pr h s [ 2 ) ) ) ( 0 ) ; 
DELA Y_LENGTH = ( int ) ( mxGet Pr ( prhs ( 3 ] ) ) [ 0 ] ; 
BG_LENGTH = ( int ) ( mxGet Pr ( prhs [ 4 ] ) ) [ 0 ) ; 
THRESHOLD = (int) (mxGetPr (prhs ( 5 ] ) ) [ 0 ] ; 
NPTS = FG_LENGTH + DELAY_L£NGTH + BG_LENGTH; 

online (ratios , fg_vals, bg_vals, x) ; 
return ; 

} 



precursorl.st 

This MATLAB script shows how one can detect a 
rhythmic discharge type precursor to seizure from 
raw data. 

% precursorl.m 

% load iir_spike_f ilt; * b and a are IIR filter coeffs 
% t~0:dt:tf; 

% iprint=0; % Set to 1 if you want print plots 

% 

% pat^'fl 1 ; eval(['load pat_' pat]); 

% x=2*(dat(: ,10)-dat(: ,ll)+dat(: ,12))-dat(: ,13)-dat(: ,9) ; 
* ^Threshold related parameters: 

% c=100; % Threshold for ratio of filter output squared to 
background 

% ncross=2 ; % Number of "rhythmic crossings'* required for 
precursor detection 

% roin_space=. 5 ; max_space=lO; % You must have at least 2 
spikes 
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% separated by at least 1 sec and less than 10 sec 
% bg_time*=5 0 ; 

% (c) Mark Frei 10-21-94 

% Remark (May 9 5) : use of sharpl.m is more effective than 
% this program for detecting these precursors. 

xlf lag«0; 

elf ;xf=f ilter (b, a, x) ;xf =xf . ~2 ; 
bg-median (xf ( 1 : 24 0*bg_time) ) ; 
tl=t (f ind(xf /bg>=c) ) ; 
dtl=diff (tl) ; 

k=f ind (min_space<=dtl & dtl<=max_space) ; 
if length (tl)>=ncross, 
k= [ l ; 1+k ( : ) ] ; 
if ( length (k) >= ncross); 

t2=tl(k(ncross: length (k) ) ) ; 
subplot (2 ,1,2); 

stem (t2 , ones (size (t2 ) ) ) ;set (gca , ' xlim ' , [ 0 
240] , 'ylim* , [0 1.5]) 

title ([ 1 Precursor detection with ' int2str ( ncross ) 

' crossings required , min_space «= * 
num2str (min_space) . . . 

• , max_space = 1 num2str (max_space) ) ) 
t_predict=t2 (1) ; 

disp ([• Precursor detection at t = ■ 
num2str (t_predict) ] ) ; 
else 

disp ('Not enough rhythmic crossings'); xlflag-l; 

end 
else 

disp ([•Only 1 int2str ( length (tl) ) * crossings of 
threshold: no detection 1 ]); return 
end 

subplot (2,1, 1) ; 
plot<t,xf/bg) 

setfgca, 'xlim' , [0 240 ] , • ylim ' , ( 0 min ( 5*c , max (xf /bg) ) ] ) 
eo-etime ( elec_onset , to ) ; co=etime (clin_onset , to) ;aeb=etime(e 
b,t0) ; 

add eventsv ( eo , co t aeb ) ; 

plot([0 240] , [c c] , 'm' ) ; 

title ([ 'Seizure prediction: pat_' pat ' => 1 
num2str (co-t__predict) . . . 

* sec prediction 1 ] ) 
if xlflag*==l, xlabel('Not enough rhythmic crossings'); 
else; xlabel( , Time (sec) ') ;end 
orient landscape; date_stamp; 
if iprint=«l, 

print -dps -P651 

end 
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precursor2 .m 

This MAT LAB script shows how one can detect the seizure 
precursor 

consisting of signal energy attenuation for a specified 

duration in 

a given input signal. 

%precursor2 .m 
% 

% This mfile gives an example of how one can detect 
quieting as 

% a precursor to seizure... 
% 

%load pat_a2; co=294 ; dt=.005; 
%n=length(dat) ; t=(0: (n-l) )*dt; tf=max(t) ; 
%nwin=f loor (tf ) -5 ; 

%chan=13 : 32 ; % Look for quieting on the grid electrodes 

(for example) 

%y=dat ( : , chan ) . * 2 ; y=sum ( y * ) 1 ; 
%c=5; % Set threshold 

% (c) Mark Frei 10-24-94 

rl=l:1000; % Look for quieting of duration 5 sec. 

Eavg=zeros (nwin, 1) ; 

% Overlapping 5 s windows translated by 1 sec each 

for j=l :nwin, Eavg( j ) =raean(y (200* ( +rl) ) ;end 

r«=l : 200*200 ; % Use 1st 200 sec to compute background 

bg=mean (y (r) ) ; % Could also use median but this weighs 

"spike quieting" heavily 

elf ;plot(4+(l:nwin) ,bg. /Eavg) 

xl=get (gca , 1 xlirc* ) ; 

hold on;plot([0 xl(2)],[c c],'m») 

t_detect=4+min(f ind (bg. /Eavg>=c) ) 

t_predict=co-t_detect 



precursor 3 -m 

This MATUVB script shows how one can detect the seizure 
precursor 

characterized by a significant drop in median frequency 
w.r.t. PSD 

for a given input signal. 

%precursor3 .m 

% 

% This mfile gives an example of how one can use an 
% abrupt drop in median frequency as a precursor to 
% seizure in some patients. 
% 



> c once, i \rt<j—z- 

BNSOOOD: <WO 9726823A1J.> 



WO 97/26823 



PCT/US97/01037 



- 145 - 

% (c) Mark Frei 10-26-94 

load pat_al % File containing data matrix dat 
C =10; ^Threshold 

dt=.005;n=length(dat) ; t= ( 0 : ( n-1 ) ) *dt ; tf=max(t); 

chan=l:32; nchan=length { chan ) ; 

nf=256; 

psd_stats 

bg=median(med(l : f loor (nwin/2 ) , : ) ) ; 
s=zeros ( length (ned) , l) ; 

for j=chan, s=s+ (med ( : , j ) /bg { j ) ) . * ( -2 ) ; end 

plot (nf *dt*(l:nwin) , (s/nchan) . *2) 

t detect=nf *dt*min (f ind ( (s/nchan) . *2>=c) ) 



Upcrossings 

This section contains the MATLAB mfiles necessary to 
compute the time of a signal's upcrossings from below a 
certain level to above that level. It can be used to 
detect signal rhythmicity and neuronal hypersynchrony if, 
e.g., the standard deviation of the most recent say 10 
inter-zero-upcrossing times is small relative to 
background. This is used to enhance seizure detection and 
artifact rejection, and can also be useful for precursor 
detection. 

upcross2 .m 

function t_upcross=upcross2 (t,x,c) 
% function t_upcross=upcross2 (t ,x, c) 
% 

% This function returns the times at which the linearly 
% interpolated function x(t) up-crosses the line~x(t)=c. 
% This version does not assume that the points of t are 
% evenly spaced. If they are, use upcrossl.m 
% 

% (c) Mark Frei 10-15-94 
i_up=f ind (dif f (x>=c)=«l) ; 

t_upcross=(c-x (i_up) ) . * (t ( i_up+l) -t ( i_up) ) . / (x(i_up+l) -x(i 
up) )+t(i_up) ; 



psd_stats 

This section contains the MATLAB mfiles necessary to 
compute statistics of the power spectral density of a 
signal such as mean, modal, and median frequency (and well 
as other \%-tiles and moments), the frequency variance, 
and average signal energy. The main routine psd\_stats.m 
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calls pdf \_ststs2 .m to compute general statistics of a 
probability density function. The "^built-in 1 ' MATLAB 
function PSD.m is used to compute the power spectral 
density {Welch's periodogram method is used here). 

%psd_stats.m 

% Need dat and t defined to run 

%chan=l : 16 ; % channels to analyze in dat 

%pat=' (pat, evnt, chan) ' ; %plot title header 
%dt=l/240; % Time step 

%nf=512; % Number of points per window (and 2* nuro 

of freqs) 

%ql= . 25 ;q2= . 7 5 ; % Used for frequency percentiles 
% (c) Mark Frei 7-20-94 

if " exist ( 'chan' ) , chan=l :min ( size (dat )) ; end 

if "exist ( 'pat* ) ,pat=' (pat , evnt , chan) * ;end 

if "exist ( 'dt* ) ,dt=l/240;end 

if "exist ( 'nf 1 ) , nf=512;end 

if "exist( 'ql' ) ,ql=.25;end 

if "exist ( 'q2 ») ,q2=. 75;end 

if "exist ( *tf 1 ) ,tf=ceil(max(t) ) ;end 

nchan* length (chan) ; 
n= length (dat) ; 

r=l:nf ;p=zeros(nf/2+l,nchan) ;nwin=f loor ( n/nf ) ; 
mu=zeros(nwin, nchan) ; mu2=mu ; sigma=mu ;med=mu ; % 
Preal location 

maxp=mu ; f _lov=mu ; f_high=mu ; Eavg=mu ; 

h_wait=waitbar (0, •Computing Statistics of PSD ...*); 
for k=l:nwin f 

rl=r+(k-l) *nf ; 

for i=l: nchan, 

[p( : , i) ,f ]=psd (dat (rl, chan (i) ) ,nf , 1/dt) ; % Welch 
periodogram method 

end 

(mu(k, : ) ,mu2 (k ( :) ,sigma(k, : ) ,med(k, : ) ,maxp(k, : ) , f low(k, : ) , 
f_high(k, :)]=... 

pdf_stats2(p,f ,ql,q2) ; 

Eavg(k, : ) =mean (dat (rl , chan) . *2) ; 

waitbar(k/nwin) ; 

end 

delete (h_wait) ; 

%pdf _stats2 . m: 
function 

( mu , mu2 , s igroa , med , raaxp , x__low , x_high ] =*pdf _stats2 ( p , x , q 1 , q2 ) ; 
%f unction 

( inu , mu2 , s igma , roed , maxp , x_low , x_high ) =pdf _stat s2 ( p , x , ql , q2 ) ; 
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e-p; 

g=0:dt: (p-l)*dt; 
f=g.'2; 
h=ones ( 1 , p) ; 
A=c*e-d~2 ; 
B=c*c-b*e ; 
C=b*d-c~2 ; 

D=a*A+2*b*c*d-e*b~ 2-c" 3 ; 
H=(A B Cj/D; 

f ilt_coeffs=2*H*(f ; g; h); 
r«p:-l:2; 

s=f ilter (f ilt_coef f s, l,y) ; 

% Now shift the output to match temporally with input 

s(l:p)=zeros(p, 1) ; 

s(i:f loor (p/2) )={ ] ; 

s= [ s ( : ) ; zeros (npts-length ( s ) , 1) ] ; 
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CLAIMS 

1 . A method of delecting the occurrence of abnormal activity in the brain oi a subject, said 
method comprising the steps of: 

(a) receiving into a signal processor input signals indicative of the subject s brain acuviiy; 

(b) determining ictal components in said input signals by applying to said input signals a first filler 
configured for extracting and enhancing icial components from said input signals; 

(c) measuring the ictal activity in a foreground epoch of said input signals by applying art order-statistic 
filter to said ictal components corresponding to said epoch to produce a foreground measure of ictal 
activity; 

(d) determining whether said foreground measure reaches a threshold level, such being indicative of 
the occurrence of said abnormal activity and 

(e) performing steps (a) - (d) while said abnormal activity is c<xurring. 

2. The method as set forth in claim 1 , step (b) including the step of using a digital filler as 

said first filter. 

3 Tne meihod as se: forth in claim 2, step (b) including the step of using z finite impulse 
response filter as said digital filter. 

4 The method as set forth in claim 2. step (b) including the step of using an infinite impulse 
response filter as said digital filter. 

5. The method as set forth in claim 1 , step (b) including the step of using an analog filter as 

said first filter. 

6. The method as set forth in claim I , step (b) including the step of squaring the results of 
the application of said first filter to produce said ictal components. 

7. The method as set forth in claim I , 

step (c) including the step of measuring the ictal activity in a background epoch of said input signals by 
applying said order- statistic to said ictal components corresponding to said background epoch to 
produce a background measure of ictal activity, said background epoch occurring before said 
foreground epoch, 

step (d) including the step of determining whether the ratio of said foreground measure to said background 
measure reaches said threshold level. 

8. The method as set forth in claim 7, step (c) including the step of configuring said 
foreground and background epochs to present a time delay therebetween. 
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% 

% Inputs: vector x and matrix p 

% such that x is increasing of length N, and 

% p is N x M, with each column of p is a 

% (possibly non-normalized) probability density 

% function p_j (x) . 

% 

% Optional: ql and g2 (low and high percentiles to use in 
returning 

k the values x_low and x_high) 

k 

* Output: 
k 

% 1 x M vectors mu, mu2 , sigma, and med, where 

% inu(j) is the mean, 

% rou2(j) is the second moment, 

% sigma(j) is the variance, and 

% med(j) is the median (50-th percentile) 

% corressponding to the j-th column of p. 

% xnaxp(j) is the x-value at which p_j (x) is maximized 

% x_low(j) and x_high(j) are defined by the relation 

k Prob(X < p_j (x_low) )=ql, Prob(X < p_j (x_high) ) =q2 

% 

k (med(j) is the x value which divides the pdf p_j(x). 

% in half in terms of area under the curve) . 

% 

% (c) Mark Frei 6-13-94 
x=x ( : ) ; 

[N,M)=size(p) ; 
if min(N,M)==l, 

p-p(: ) ; 

else 

if N "= length (x) , 

disp(' Error, number of rows in p should .equal the 
length of x ' ) ; 
end 

nu« zeros (1,M) ;rau2=mu; sigroa«mu ;med=mu; 
if nargout > 4 , 

maxp=mu ; x_low=mu ,*x_high=mu ; %preal locate 

end 

end 

k Convert to pdf(s) 
s=trapz(x,p) ; 
for j=f ind (s" = l ) ; 

PC.j)-pC,j)/s(j) ; 
end 

% Compute stats 
x2=x . * 2 ; 
for j=l:M, 

*u(j)=trapz(x,x.*p(: , j)) ; 

mu2 (j)=trapz(x,x2.*p(: , j ) ) ; 

end 
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sigma=rau2-mu . " 2 ; 

if nargout > 4 , 
(a,b)=max(p) ; 
maxp=x(b) ' ; 

end 

. dx=diff (x) ;k«l: (N-l) ; 
for j=l:M, 

ipdx=[0; . 5*cumsuro( (p (k, j ) +p (k+1 , j ) ) . *dx) ] ; 

ind=[0; f ind(dif f ( ipdx) >0) )+l ; %makes ipdx increasing 

med( j)=tablel([ipdx(ind) x(ind) ] , .5) ; 

if nargin > 2 , 

x_low(j)=tablel([ipdx(ind) x(ind) ] ,ql) ; 

end 

if nargin > 3 , 

x_high( j)=tablel( (ipdx(ind) x(ind) ] ,q2) ; 

end 

end 
end 

sharp 1 -m 

function s=sharpl (y , p , dt) ; 
% function s=sharpl (y , p, dt) ; 
% 

k This function returns the 2nd-der ivative 

% (i.e. acceleration) of the best approximating 

% parabola (LS sense) to data y, in a moving 

% window of p points. 

% 

% Warning: The FIR filter used here instills 

% a delay of approximately dt*(p/2) in time 

% (one needs approx. p/2 future points to 

% evaluate the acceleration/curvature at the 

% current time point) - We shift the output 

% backward in time to correct for this, so 

% ' that the output "sharpness," s, is the same 

% length as the input, and corresponds temporally 

% with where the sharpness occured (not when it 

% was available from computations (which would 

* be about (p/2)*dt later...)- 

H (c) Mark G. Frei 6-11-95 

y=y(0 ; 

npts=length(y) ; 

if nargin <3 , dt=l ; end 

if nargin <2 , p=3 ; dt=l ; end 

a=(p*(p-l)*(2*p-l)*(3*p"2-3*p-l) *dt~4) /30; 
b=(dt"3*(p-l) -2*p~2) /4 ; 
c=(p*(p-l)*(2*p-l)*dt-2) /6; 
d=p*(p-l)*dt/2; 
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9 The method as set forth in claim 8, step (c I including the step confipunng said foreground 
epoch as about two seconds, said background epoch as about twenty seconds and said delay as about one second 

10 The method as set forth in claim 7. step (c) including the step of continuously updating 
said foreground and background measures. 

1 1 The method as set forth in claim 10, said threshold level being a first threshold level, step 

(c) including the step of suspending updating said background measure if said ratio reaches a second threshold level. 

12. The method as set forth in claim 1 , step (e) including the step of performing steps (a) - 

(d) before the onset of the clinical component of the seizure thereby being predictive of the clinical component 

13. The method as set forth in claim 1 2, step (e) including the step of performing steps (a) - 
(d) before the onset of the eieclrographic component of the seizure thereby being predictive of the electrographic 
component 

14 The method as set forth in claim I , step (b) including the step of selecting said first filter 
from a filter bank including o plurality of filters configured for extracting and enhancing said ictal components. 

1 5. The method as set forth in claim 14, step (b) including the step of selecting said first filter 
as the filter from said filter bank providing the greatest differentiation of ictal components. 

16. The method as set forth in claim 1 5, said signal processor including memory means for 
storing said filter bank, step (b) including the step of retrieving a plural iry of filters from said filter bank and applying 
said plurality of filters to said input signals in said signal processor and therein selecting said first filter as the filter 
from said plural iry providing the greatest differentiation of ictat components. . 

1 7. The method as set forth in claim 1 , further including the step of configuring said first filter 
in order to increase the differentiation of ictal components for the subject. 

1 8. The method as set forth in claim 1 , 

stcp(b) including ihe steps of using one of a finite impulse response filter and an infinite impulse response 
filter as said digital filter, and squaring the results of the application of said first filler to enhance 
said ictal components, 

step (c) including the steps of measuring the ictal activity in a background epoch of said input signals by 
applying said order-statistic to said ictal components corresponding to said background epoch to 
produce a background measure of ictal activity and continuously updating said foreground and 
background measures, said background epoch occurring before said foreground epoch and said 
foreground and background epochs presenting a time delay therebetween. 
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step (d> including ihe step of determining whether the ratio ol said foreground measure 10 said background 

measure reaches said threshold level, 
said threshold level being a first threshold level, step (c) including the step ol" suspending updating said 
background measure tf said ratio reaches a second threshold level. 

19. The method as set fonh in claim I , step (a) including the step of receiving said signals 
from at least one electrode operable for detecting the subject's brain activity and for producing said signals indicative 
thereof, said at least one electrode being selected from the group consisting of a scalp electrode and an implanted 
electrode. 

20. The method as set forth in claim 1 , step (a) including the step of receiving said signals 
from a memory device. 

2 1 The method as set forth in claim I . step (b) including the step of analyzing said signals 
1 5 in said signal processor selected from the group consisting of a microprocessor and a computer. 

22. The method as set forth in claim 1 , step (b) including the step of analyzing said signals 
in a signal processor implanted within the subject. 

20 23. The method as set fonh in claim 1 , step <b) including the steps of analyzing said signals 

by using a wavelet filter as said first filler to determine corresponding wavelet coefficients of said signals, using said 
wavelet coefficients to determine a power density distribution, and comparing said power density distribution with 
said threshold level wherein the crossing of said threshold leve) by said distribution indicates the occurrence of the 
seizure. 

25 

24 The method as set forth in claim I . step (b) including the step of analyzing said signals 
using windowed Fourier and inverse Fourier transforms as said first filter. 

25. The method as set forth in claim 1 including the step of producing an output in response 
30 lo the indication of the occurrence of a seizure as said abnormal activity with said output taken from the group 

consisting of administering a medicament, electrically stimulating a portion of the subject's brain, magnetically 
stimulating a portion of the subject's brain, inhibiting activity in a portion of the subject *s brain, electrically stimulating 
a nerve of the subject, recording said signals, activating an alert, stimulating physiological receptors of the patient, 
heating at least a portion of the subject's brain, cooling at least a portion of the subject's brain, facilitating activity in 
35 a portion of a subject's brain, disfacilitating activity m a portion of a subject's brain, and ablating a portion of the 
subject's brain. 

26. The method as set forth in claim 25. further including the step of using a device implanted 
with the subject for producing said output. 

40 
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27. The method as set fonh in claim I , said first filler including a genetically designed filler 

28 A method of predicting the occurrence of a seizure in ihe brain of a subject, said method 
comprising the steps of 

(a) receiving into a signal processor input signals indicative of the subject s brain activity; 

(b) analyzing said signals for at least one precursor predictive of the occurrence of a seizure in the 
subject; and 

(c) upon occurrence of said at least one precursor, producing an output in response before the 
occurrence of the seizure. 



29. The method as set forth in claim 28, step (b) including the step of delecting epileptiform 
discharges in said signals, determining the sharpness of each of said discharges by determining a parabola of optimal 
fit for each of said spikes, determining an index of relative sharpness of each of said discharges by comparing each 
sharpness to the sharpness of other discharges in a time window, determining whether said indc\ reaches said 

1 5 predetermined level, said spike presenting a polarity, and determining whether said spikes reaching said predetermined 
level fit a predetermined pattern, such being predictive of a seizure. 

30. The method as set forth in claim 29, step (b) including the step of determining said 
precursor by detecting epileptiform spikes by using a signal analysis filter to extract spike shape coefficients, 

20 determining a ratio of spike shape coefficients squared lo background signal informal ion, determining whether said 
ratio exceeds said predetermined level and if so, determining whether spikes exceeding said ratio fit a pattern 
determined as being predictive of a seizure 

3 1 The method as set fonh in claim 28, step tb) including the step of determining said 
25 precursor by determining a ratio of current signal energy compared to background energy and determining whether 
said ratio exceeds said predetermined level, such being predictive of a seizure. 

32. The method as set forth in claim 28, step (b) including the step of determining said 
precursor by determining a ratio of median frequency to background median frequency and determining whether said 
30 ratio exceeds said predetermined level, such being predictive of a seizure. 

33 A method of analyzing the activity slate of the brain of a subject, said method comprising 

the steps of: 

(a) receiving into a signal processor signals representative of the subject s brain activity. 
35 (b) analyzing said signals in said signal processor for presence of information indicating the current 

activity state of the subject's brain; 

(c) accomplishing step (b) while said activity state is occurring; and 

(d) producing an output in response to said presence in said information 
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34. The method as sei forth in claim 1 . step (b) including the step of using a filter as said first 
filter taken from the group consisting of a digital filter, a nonlinear filter, and adaptive filler, a correlation integral, an 
arc length differential and a temperature filter. 

5 35. The method as set forth in claim 1 , said method including the method of detecting the 

occurrence of an epileptic seizure as said seizure. 

36. The method as set forth in claim I further including the step of receiving other biological 
signals concerning the subject and using said biological signals in detecting the occurrence of a seizure, said biological 

1 0 signals being representative of biological functions of the subject selected from the group consisting of respiratory 
activity, concentrations of glucose, free radicals, neurotransmitters, blood, body tissues, brain temperature, heart 
activity, muscle activity, ocular aciiviry. magnetic fields, skin resistance, and electrical fields. 

37. The method as set forth in claim I further including the step of providing biofeedback to 
1 5 the subject indicative of a seizure or epileptiform discharge. 

38. The method as set forth in claim 1 , including the step of measuring said ictal activity in 
said foreground epoch against a background signal generated using time and slate -weighted averaging. 

20 39. The method as set forth in claim 38, including the step of using an exponentially forgetting 

tune averaging as said time and state-weighted averaging. 

40 A seizure detection, prediction and treatment apparatus for detecting the occurrence of 
a seizure in the brain of a subject, said method comprising: 

receiving means for receiving input signals indicative of the subject's brain activity; 
signal processing means coupled with said receiving means for 

determining ictal components in said input signals by applying to said input signals a first filter 

configured for extracting and enhancing ictal components from said input signals, 
measuring the ictal acuviry in a foreground epoch of said input signals by applying an order- 
statistic filter to said ictal components corresponding to said epoch to produce a 
foreground measure of ictal activity, and 
determining whether said foreground measure exceeds a predetermined level, such being indicative 

of the occurrence of a seizure, and 
output means for producing an output indicating the occurrence of the seizure while the seizure is 
occurring. 



30 



BNSDOCID. <WO S726623A1_I, > 



WO 97/26823 



PCT7US97/01037 



1/15 




+ 



SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



2/15 



PCT/US9 7/0103 7 



Fig. 2 



CURRENT 
TIME 

i 

— 20 SEC. >|- -1 SEC . 2 SEC. > 

BACKGROUND j DELAY ] FOREGROUND] TIME 

WINDOW i WINDOW i WINDOW i 



c-» idctiT! frzr c t_'rrc"r / o » » « rr pc^ 

BNSDOClD: <WO 9726S23A1. I.> 



WO 97/26823 



PCT/US97/01037 



Fig.3A GENERIC FIR FILTER 

0.6 I 1 > ' — -i 1 1 1 




_q 4 I 1 i i i i i i 1 i 

" 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

TIME (SEC.) 



Fig.3B 

0.01 5 | 1 1 1 1 1 1 1 r 



0.01 - 



0.005 - 



0 




0 5 10 15 20 25 30 35 40 45 50 

FREQUENCY (HZ.) 



SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



PCTAJS97/01037 



Fig. 4 



4/15 



ECoG CONTAINING A SEIZURE EVENT 




-0.3 



0 



60 



70 



80 



90 



100 




> 
E 



> 
E 




180 



190 



200 210 220 

TIME (SEC.) 



230 



240 



SUBSTITUTE SHEET (RULE 26> 



SNSDOCIO: <WO..„9726823A1.l > 



WO 97/26823 



PCT/US97/01037 



5/15 




SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 PCT/US97/01 037 

6/15 




BNSDOClD: <WO S726823A1 J .> 



SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



PCT/US97/01037 




SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



PCT/US97/01037 




bNSDOCID. <WO &726823A1. I > 



WO 97/26823 



PCT/US97/01037 



9/15 

_ Fig. 9 

PSD OF MR-FILTER FOR SPIKE PRECURSOR 

0.3 I 1 1 r 1 1 1 1 1 1 1 




0 2 4 6 8 10 12 14 16 18 20 



FREQUENCY (HZ.) 



SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



PCTYUS97/01037 



10/15 



ECoG DATA SHOWING A PRECURSOR 




-0.4 



80 82 84 86 88 90 92 94 96 98 100 




-0.4 



100 102 1 04 1 06 108 110 112 114 116 118 120 



> 
E 




-0.4 



120 122 124 126 128 130 132 134 136 138 140 



0.4 



> 
E 



-0.4 



1 1 1 1 1 T 


i< — CLINICAL 




i ONSET 


1 L 1 1 I 


.!_ - 1 J 



140 142 144 146 148 150 152 154 156 158 160 

TIME (SEC.) 



C» IDCTfT! 1 



BNSDOCID <WO . 9726823A1. 1. > 



WO 97/26823 



PCT/US97/01037 



11/15 




SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 




I 



SUBSTITUTE SHEET (RULE 26) 



BNSDOCID <WO„ 9726823A1 I. > 



WO 97/26823 



PCT/US97/01037 




SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



Fig. 14A 



Fig. 14B 



0.5- 



0.4- 



wo.a- 

0l 



02- 



0.1 




14/15 



PCT/US97/01037 



1024 POINTS OF ECoG SIGNAL 




2 2.5 
TIME (SEC.) 



- -MODAL FREQUENCY 

- - MEDIAN FREQUENCY 
MEAN FREQUENCY 




4 5 6 7 
FREQUENCY (HZ.) 



8 



10 



BNSDOCJD: <WO 9726823A1_1_> 



SUBSTITUTE SHEET (RULE 26) 



WO 97/26823 



PCT/US97/01037 



Fig. 15A 

0.05 



15/15 

ECoG DATA 




Fig. 1 5B ABSOLUTE VALUE OF WAVELET COEFFICIENTS 



l! i. nil I ill 



' i .ii-.i.ii.i. . i.h . iiiii., i. i. H.t i . i in M it. . .i.ii.ii.i.ii 



In Hlh, | 



LU 

> 3 

1X1 



I I l I 



i . 1 1 



0 0.2 0.4 0.6 



0.8 1 1.2 
TIME (SEC.) 



1.4 1.6 1.8 



SUBSTITUTE SHEET (RULE 26) 



INTERNATIONAL SEARCH REPORT 



International application No. 
PCT/US97/01037 



A. CLASSIFICATION OF SUBJECT MATTER 

IPC(6) :A61B 5/0476; A61N 1/18, 36 
US CL .128/731, 732; 607/45 
According to International Patent Classification (IPC) or to both national classification and IPC 


B. FIELDS SEARCHED 






Minimum documentation searched (classification system followed by classification symbols) 
U.S. : 128/731. 732; 607/45 


Documentation searched other than minimum documentation to the extent that such documents arc included in the fields searched 


Electronic data base consulted during the interna tiona] search (name of data base and. where practicable, search terms used) 
APS 


C. DOCUMENTS CONSIDERED TO BE RELEVANT 


Category* 


Citation of document, with indication, where appropriate, of the relevant passages 


Relevant to claim No. 


X 


US 5,361,773 A (IVES) 08 November 1994, Abstract; col. 
1, lines 26-53; col. 6, lines 1-61; and claims 1-26. 


1-40 


X 


US 5,215,086 A (TERRY, JR. et at) 01 June 1993, Abstract; 
col. 1, lines 29-47; col. 11, lines 49 + ; and claim 1-18. 


1-40 


Y, P 


US 5,517,251 A (RECTOR et al) 14 May 1996, col. 2, lines 
9-11; and col. 6, lines -15-25. 


1-40 


| | Further documents arc listed in the continuation of Box C. [ | See patent family annex. 


" Special categories of eked documents: 

*A* document defuunf the general state of (heart which is notooMidered 
to be of particular ickvuct 

'E* earlier document published on or after the international filing date 

*L* document which may throw doubt* oo priority caumfs) or which is 
cited lo establish the pubbcanoo date of another cnaboo or other 
•nccial reason (a* specified) 

*0* document refcrrm*. to an oral disclosure, use. exhibiboa or other 
mean* 

"1*" document publiahed prior to the aUcrnauonal filing date but aster than 
the priority date claimed 


*T* later document published after the international filing date or priority 
dale and Dot m conflict with die application but cited to understand the 
praacaple or theory underlying the invention 

*X* document of paruodar relevance: the claimed invention cannot be 
cooaidered novel or cannot be considered to involve an inventive atep 
when the document ■ uken alone 

"Y" document of particular relevance; the claimed invention cannot be 
considered lo involve an inventive atcp when the document ■* 
combined with one or more other aucb document!, audi combination 
being obvious to a person skilled in the an 

'Sl' document memher of the same patent family 


Date of the actual completion of the international search 
10 APRIL 1997 


Date of mailing of the international search report 

0 7^ AY 1997 


Name and mailing address of the ISA/US 
Commissioner of Patent* and Trademarks 
Box PCT 

Washingmn. D.C. 20231 — 
Facsimile No. (703) 305-3230 ( 


^/STEPHEN HUANG 

Telephone No. (703) 308-3399 



